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Abstract 

Rotating relativistic stars have been studied extensively in recent years, both 
theoretically and observationally, because of the information one could obtain 
about the equation of state of matter at extremely high densities and because they 
are considered to be promising sources of gravitational waves. The latest the- 
oretical understanding of rotating stars in relativity is reviewed in this updated 
article. The sections on the equilibrium properties and on the nonaxisymmetric 
instabilities in /-modes and r-modes have been updated and several new sections 
have been added on analytic solutions for the exterior spacetime, rotating stars in 
LMXBs, rotating strange stars, and on rotating stars in numerical relativity. 
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1 Introduction 



Rotating relativistic stars are of fundamental interest in physics. Their bulk proper- 
ties constrain the proposed equations of state for densities larger than nuclear density. 
Accreted matter in their gravitational field undergoes high-frequency oscillations that 
could become a sensitive probe for general relativistic effects. Temporal changes in the 
rotational period of millisecond pulsars can also reveal a wealth of information about 
important physical processes inside the stars or of cosmological relevance. In addition, 
rotational instabiUties can produce gravitational waves, the detection of which would 
initiate a new field of observational asteroseismology of relativistic stars. 

There exist several independent numerical codes for obtaining accurate models of 
rotating neutron stars in full general relativity, including one that is freely available. 
One recent code achieves near machine accuracy even for uniform density models near 
the mass-shedding limit. The uncertainty in the high-density equation of state still 
allows numerically constructed maximum mass models to differ by as much as a factor 
of two in mass, radius and angular velocity, and a factor of eight in the moment of 
inertia. Given these uncertainties, an absolute upper limit on the rotation of relativistic 
stars can be obtained by imposing causality as the only requirement on the equation of 
state. It then follows that gravitationaUy bound stars cannot rotate faster than 0.28 ms. 

In rotating stars, nonaxisymmetric perturbations have been studied in the Newto- 
nian and post-Newtonian approximations, in the slow-rotation limit and in the Cowl- 
ing approximation but fully relativistic quasi-normal modes (except for neutral modes) 
have yet to be obtained. A new method for obtaining such frequencies is the time- 
evolution of the full set of nonhnear equations and frequencies of quasi-radial modes 
have already been obtained. Time-evolutions of the linearized equations have also im- 
proved our understanding of the spectrum of axial and hybrid modes in relativistic 
stars. 

Nonaxisymmetric instabilities in rotating stars can be driven by the emission of 
gravitational waves (CFS-instability) or by viscosity. Relativity strengthens the former, 
while it weakens the latter. Nascent neutron stars can be subject to the / = 2 bar mode 
CFS-instability, which would turn them into a strong gravitational wave source. 

Axial fluid modes in rotating stars (r-modes) received considerable attention since 
it was discovered that they are genericaUy unstable to the emission of gravitational 
waves. The r-mode instabiUty could slow down newly-born relativistic stars and limit 
their spin during accretion-induced spin-up, which would explain the absence of mil- 
lisecond pulsars with rotational periods less than ~ 1 .5 ms. Gravitational waves from 
the r-mode instability could become detectable if the amplitude of r-modes is of order 
unity. Recent 3D simulations show that this is possible on dynamical timescales, but 
nonhnear effects seem to set a much smaller saturation amphtude on longer timescales. 
Still, if the signal persists for a long time (as has been found to be the case for strange 
stars) even a small amplitude could become detectable. 

Recent advances in numerical relativity have enabled the long-term dynamical evo- 
lution of rotating stars and several interesting phenomena, such as dynamical instabili- 
ties, pulsations modes, neutron star and black hole formation in rotating collapse have 
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now been studied in full general relativity. The current studies are limited to relativistic 
polytropes, but new 3D simulations with reaUstic equations of state should be expected 
in the near future. 

The present article aims at presenting a summary of theoretical and numerical meth- 
ods that are used to describe the equilibrium properties of rotating relativistic stars, 
their oscillations and their dynamical evolution. It focuses on the most recently avail- 
able preprints, in order to rapidly communicate new methods and results. At the end of 
some sections, the reader is pointed to papers that could not be presented in detail here 
or to other review articles. As new developments in the field occur, updated versions 
of this article will appear. 

2 The Equilibrium Structure of Rotating Relativistic 
Stars 

2.1 Assumptions 

A relativistic star can have a comphcated structure (sohd crust, magnetic field, possible 
superfluid interior, possible quark core etc.). Still, its bulk properties can be computed 
with reasonable accuracy by making several simplifying assumptions. 

The matter can be modeled to be a perfect fluid because observations of pulsar 
gUtches have shown that the departures from a perfect fluid equiUbrium (due to the 
presence of a sohd crust) are of order 10^^ (see [109]). The temperature of a cold 
neutron star does not affect its bulk properties and can be assumed to be K, because 
its thermal energy (<< IMeV 10^" K) is much smaUer than Fermi energies of the 
interior (> 60 MeV). One can then use a zero-temperature, bawtropic equation of state 
(EOS) to describe the matter: 

e = e(P), (1) 

where 8 is the energy density and P is the pressure. At birth, a neutron star is expected 
to be rotating differentially, but as the neutron star cools, several mechanisms can act 
to enforce uniform rotation. Kinematical shear viscosity is acting against differential 
rotation on a timescale that has been estimated to be [98, 99, 76J 

where p, T and R are the central density, temperature and radius of the star. It has also 
been suggested that convective and turbulent motions may enforce uniform rotation on 
a timescale of the order of days [150]. In recent work, Shapiro [268] suggests that 
magnetic braking of differential rotation by Alfven waves could be the most effective 
damping mechanism, acting on short timescales of the order of minutes. 

Within roughly a year after its formation, the temperature of a neutron star becomes 
less than 10^ K and its outer core is expected to become superfluid (see [225] and 
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references therein). Rotation causes superfluid neutrons to form an array of quantized 
vortices, with an intervortex spacing of 



(3) 



where ^2 is the angular velocity of the star in 10^ s^ ^ On scales much larger than the 
intervortex spacing, e.g. on the order of 1 cm, the fluid motions can be averaged and 
the rotation can be considered to be uniform [286]. With such an assumption, the error 
in computing the metric is of order 



assuming R ^ 10 km to be a typical neutron star radius. 

The above arguments show that the bulk properties of an isolated rotating relativis- 
tic star can be modeled accurately by a uniformly rotating, zero-temperature perfect 
fluid. Effects of differential rotation and of finite temperature need only be considered 
during the first year (or less) after the formation of a relativistic star. 

2.2 Geometry of Spacedme 

In general relativity, the spacetime geometry of a rotating star in equilibrium can be 
described by a stationary and axisymmetric metric gab of the form 



where V, (0 and /u are four metric functions which depend on the coordinates r 
and 9 only (see e.g. Bardeen & Wagoner [26]). Unless otherwise noted, we will 
assume c = G = 1 . In the exterior vacuum, it is possible to reduce the number of metric 
functions to three, but as long as one is interested in describing the whole spacetime 
(including the source-region of nonzero pressure), four different metric functions are 
required. It is convenient to write in the the form 



where B is again a function of r and 9 only [24] . 

One arrives at the above form of the metric assuming that i) the spacetime has a 
timelike Killing vector field t" and a second Killing vector field (j)" corresponding to 
axial symmetry, ii) the spacetime is asymptotically flat, i.e. taf^ = — 1, '^a'^" = +°° ™d 
ta<^" = at spatial infinity. According to a theorem by Carter [21], the two Killing 
vectors commute and one can choose coordinates x" = ? and = (|) (where x", a = 0. .3 
are the coordinates of the spacetime), such that t" and are coordinate vector fields. If, 
further, the source of the gravitational field satisfies the circularity condition (absence 
of meridional convective currents), then another theorem [56] shows that the 2-surfaces 
orthogonal to t" and can be described by the remaining two coordinates x^ and x^. 




(4) 



ds^ = -e^^dt^ + _ (Mtf + e^f{dr^ + r^dQ^), 



(5) 



(6) 
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A common choice for and are quasi-isotropic coordinates, for which gr^ = 0, 
gOQ = r^grr (in spherical polar coordinates), or g^^ = 0) = r^gtm (in cylindrical 
coordinates). In the slow-rotation formalism by Hartle [140], a different form of the 
metric is used, requiring ^ee = gi/ii / sin^ 6. 

The three metric functions V, \|/ and CO can be written as invariant combinations of 
the two KilUng vectors t" and (j)", through the relations 

tj" = git, (7) 

^a<f = g^, (8) 
tar = gt^, (9) 

while the fourth metric function p determines the conformal factor e^^' that character- 
izes the geometry of the orthogonal 2-surfaces. 

There are two main effects that distinguish a rotating relativistic star from its non- 
rotating counterpart: the shape of the star is flattened by centrifugal forces (an effect 
that first appears at second order in the rotation rate) and the local inertial frames are 
dragged by the rotation of the source of the gravitational field. While the former effect 
is also present in the Newtonian limit, the latter is a purely relativistic effect. The study 
of the dragging of inertial frames in the spacetime of a rotating star is assisted by the in- 
troduction of the local Zero- Angular-Momentum-Observers (ZAMO) [23, 24]. These 
are observers whose worldlines are normal to the t = const, hypersurfaces and are also 
called Eulerian observers. Then, the metric function (O is the angular velocity of the 
local ZAMO with respect to an observer at rest at infinity. Also, is the time dilation 
factor between the proper time of the local ZAMO and coordinate time t (proper time at 
infinity), along a radial coordinate line. The metric function \|/ has a geometrical mean- 
ing: is the proper circumferential radius of a circle around the axis of symmetry. In 
the non-rotating limit, the metric (5) reduces to the metric of a non-rotating relativistic 
star in isotropic coordinates (see [323] for the definition of these coordinates). 

In rapidly rotating models, an ergosphere can appear, where gn > 0. In this region, 
the rotational frame-dragging is strong enough to prohibit counter-rotating time-like 
or null geodesies to exist and particles can have negative energy with respect to a sta- 
tionary observer at infinity. Radiation fields (scalar, electromagnetic or gravitational 
waves) can become unstable in the ergosphere [104], but the associated growth time is 
comparable to the age of the universe [66]. 

The asymptotic behaviour of the metric functions V and CO is 

V ~ -f cose , (10) 

r r^ 

2J 

(0 ~ -, (11) 

r^ 

where M, J and Q are the gravitational mass, angular momentum and quadrupole mo- 
ment of the source of the gravitational field (see Section 2.5 for definitions). The 
asymptotic expansion of the dragging potential CO shows that it decays rapidly far from 
the star, so that its effect will be significant mainly in the vicinity of the star. 
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2.3 The Rotating Fluid 



When sources of non-isotropic stresses (such as a magnetic field or a sohd state of 
parts of the star), viscous stresses and heat transport are neglected in constructing an 
equilibrium model of a relativistic star, then the matter can be modeled as a perfect 
fluid, described by the stress-energy tensor 

T'''' = {E + P)u''u''+Pg'''', (12) 

where u" is the fluid's 4-velocity. In terms of the two KilUng vectors t" and (j)", the 
4- velocity can be written as 

u"=^l=r+nr), (13) 

V 1 — 

where v is the 3-velocity of the fluid with respect to a local ZAMO, given by 

v=(r2-to)eV-\ (14) 

and Q. = u'^/u' = d(^/dt is the angular velocity of the fluid in the coordinate frame, 
which is equivalent to the angular velocity of the fluid as seen by an observer at rest at 
infinity. Stationary configurations can be differentially rotating, while uniform rotation 
{Q, = const.) is a special case (see Section 2.5). 



2.4 Equations of Structure 

Having specified an equation of state of the form e = e(P), the structure of the star is 
determined by solving four components of Einstein's gravitational field equations 



Rab = 87t ( Tab - \sabT 



(15) 



(where Rab is the Ricci tensor and T = Tg") and the equation of hydrostationary equi- 
librium. Setting t^=fi+v, one conmion choice for the gravitational field equations is 
[54] 



V-(BVv) 
V- (rsineVB) 



ir^sin^eS^e-'^^VoVo) 
2 



+47tfie2C-2v 



(e + P)(l+v2) 



l-v2 

^2,2C-4v (£ + -P)v 
1-V2 ' 

167trsineBe^^-2^P, 



2P 



= -167trsineB^e^ 



(16) 

(17) 

(18) 
(19) 



supplemented by a first-order differential equation for ^ (see [54]). Above, V is the 3- 
dimensional derivative operator in a flat 3-space with spherical polar coordinates r, 6, (|). 
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Thus, three of the four gravitational field equations are elliptic, while the fourth 
equation is a first-order partial-differential equation, relating only metric functions. 
The remaining non-zero components of the gravitational field equations yield two more 
elliptic equations and one first-order partial-differential equation, which are consistent 
with the above set of four equations. 

The equation of hydrostationary equilibrium follows from the projection of the 
conservation of the stress-energy tensor normal to the 4-velocity, {S^b + u'^Uh)VaT^^ = 
0, and is written as 



Pn + {£ + P) 



1 ^ ^2 
-VV,i+V 



l-v2V " Q.-0) 



0, (20) 



where a comma denotes partial differentiation and i=l.3. When the equation of state 
is barotropic then the hydrostationary equilibrium equation has a first integral of motion 

fP dP 

/ \n{u"Vat)+ / F{D.)d£l = constant = v\r^ie, (21) 

Jo t + P J sic 

where F{£1) = Uisfu' is some specifiable function of £1 only and fic is the angular ve- 
locity on the symmetry axis. In the Newtonian limit, the assumption of a barotropic 
equation of state impUes that the differential rotation is necessarily constant on cylin- 
ders and the existence of the integral of motion (21) is a direct consequence of the 
Poincare-Wavre theorem (which implies that when the rotation is constant on cylin- 
ders, the effective gravity can be derived from a potential, see [303]). 



2.5 Rotation Law and Equilibrium Quantities 

A special case of rotation law is uniform rotation (uniform angular velocity in the 
coordinate frame), which minimizes the total mass-energy of a configuration for a given 
baryon number and total angular momentum [48, 144]. In this case, the term involving 
F{£1) in (21) vanishes. 

More generally, a simple choice of a differential-rotation law is 

^/^x .2/^ (n-co)r2sin2ee2(P-v) 
^ ' ^ " ' l-(a-(0)Vsin2ee2(p-v)' 

where A is a constant [181, 182]. When A ^ 00, the above rotation law reduces to 
the uniform rotation case. In the Newtonian limit and when A — > 0, the rotation law 
becomes a so-called 7— constant rotation law (specific angular momentum constant 
in space), which satisfies the Rayleigh criterion for local dynamical stabihty against 
axisymmetric disturbances (j should not decrease outwards, dj/dO. < 0). The same 
criterion is also satisfied in the relativistic case [182]. It should be noted that differ- 
entially rotating stars may also be subject to a shear instability that tends to suppress 
differential rotation [337]. 
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Table 1: Equilibrium properties 



circumferential radius R = e^ 
gravitational mass M = f {Tab - 1 /'^■8ab)Tt''n''dV 

baryon mass Mq = / puaft^dV 

internal energy U = J uuaffdV 

proper mass Mp =Mq + U 

gravitational binding energy W = M — Mp — T 

angular momentum J = J Tab<^''n''dV 

moment of inertia I = J/D, 

kinetic energy T = 1 /2JQ. 



The above rotation law is a simple choice that has proven to be computationally 
convenient. More physically plausible choices must be obtained through numerical 
simulations of the formation of relativistic stars. 

Equilibrium quantities for rotating stars, such as gravitational mass, baryon mass, 
angular momentum etc. can be obtained as integrals over the source of the gravitational 
field. A list of the most important equilibrium quantities that can be computed for 
axisymmettic models, along with the equations that define them, is displayed in Table 
1. There, p is the rest-mass density, u = e— pc^ is the internal energy density, n" = 
^at /\^ht^^t\^^^ is the unit normal vector field to the t = const, spacelike hypersurfaces 
and dV = \/\^g\ d^x is the proper 3-volume element (with being the determinant 
of the 3-metric). It should be noted that the moment of inertia cannot be computed 
directly as an integral quantity over the source of the gravitational field. In addition, 
there exists no unique generalization of the Newtonian definition of the moment of 
inertia in general relativity and I = J/Q.isa connmon choice. 

2.6 Equations of State 
2.6.1 Relativistic Polytropes 

An analytic equation of state that is commonly used to model relativistic stars is the 
adiabatic, relativistic polytropic EOS of Tooper [313] 

P = Kp^, (23) 

e = pc2 + ^, (24) 

where K and F are the polytropic constant and polytropic exponent, respectively. No- 
tice that the above definition is different from the form P = Ke^ (also due to Tooper 
[312]) that has also been used as a generalization of the Newtonian polyttopic EOS. 
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Instead of F, one often uses the polytropic index A^, defined through 



r = i + i. (25) 

For the above equation of state, the quantity c^^^^^^^^^^^^/kV(J^~W/G has units of 
length. In gravitational units (c = G = 1), one can thus use as a fundamental 
length scale to define dimensionless quantities. Equilibrium models are then character- 
ized by the polytropic index and their dimensionless central energy density. Equilib- 
rium properties can be scaled to different dimensional values, using appropriate values 
for K. FovN <l.O(N>l .0) one obtains stiff (soft) models, while for A? - 0.5 - 1 .0, 
one obtains models with bulk properties that are comparable to those of observed neu- 
tron star radii and masses. 

Notice that for the above polytropic EOS, the polytropic index F coincides with the 
adiabatic index of a relativistic isentropic fluid 

_e + PdP 

r = Fad = ^^. (26) 

This is not the case for the polytropic equation of state P = Ke^, which satisfies (26) 
only in the Newtonian limit. 



2.6.2 Hadronic Equations of State 

The true equation of state that describes the interior of compact stars is, still, largely 
unknown. This comes as a consequence of our inabiUty to verify experimentally the 
different theories that describe the strong interactions between baryons and the many- 
body theories of dense matter, at densities larger than about twice the nuclear density 
(i.e. at densities larger than about 5 x lO^'^gr/cm^). 

Many different so-called realistic EOSs have been proposed to date which all pro- 
duce neutron star models that satisfy the currently available observational constraints. 
The two most accurate constraints are that the EOS must admit nonrotating neutron 
stars with gravitational mass of at least 1 .44M© and allow rotational periods at least 
as small as 1.56 ms, see [242, 184]. Recently, the first direct determination of the 
gravitational redshift of spectral Unes produced in the neutron star photosphere has 
been obtained [72]. This determination (in the case of the low-mass X-ray binary EXO 
0748-676) yielded a redshift of z = 0.35 at the surface of the neutron star, correspond- 
ing to a mass to radius ratio of M/R = 0.23 (in gravitational units), which is compatible 
with most normal nuclear matter EOSs and incompatible with some exotic matter EOS. 

The theoretically proposed EOSs are qualitatively and quantitatively very different 
from each other. Some are based on relativistic many-body theories while others use 
nonrelativistic theories with baryon-baryon interaction potentials. A classic collection 
of early proposed EOSs was compiled by Arnett and Bowers [18], while recent EOSs 
are used in Salgado et al. [262] and in [81]. A review of many modem EOSs can be 
found in a recent article by Haensel [135]. Detailed descriptions and tables of several 
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modern EOSs, especially EOSs with phase transitions, can be found in Glendenning's 
book [122]. 

High density equations of state with pion condensation have been proposed by 
Migdal [226] and Sawyer and Scalapino [265]. The possibility of Kaon condensa- 
tion is discussed by Brown and Bethe [50] (but see also Pandharipande et al. [240]). 
Properties of rotating Skyrmion stars have been computed in [236]. 

The realistic EOSs are supplied in the form of an energy density vs. pressure table 
and intermediate values are interpolated. This results in some loss of accuracy because 
the usual interpolation methods do not preserve thermodynamical consistency. Swesty 
[302] devised a cubic Hermite interpolation scheme that does preserve thermodynam- 
ical consistency and the scheme has been shown to indeed produce higher accuracy 
neutron star models in Nozawa et al. [235]. 

Usually, the interior of compact stars is modeled as a one-component ideal fluid. 
When neutron stars cool below the superfluid transition temperature, the part of the star 
that becomes superfluid can be described by a two-fluid model and new effects arise. 
Andersson and Comer [7] have recently used such a description in a detailed study of 
slowly rotating superfluid neutron stars in general relativity, while first rapidly rotating 
models are presented in [247]. 

2.6.3 Strange Quark Equations of State 

Strange quark stars are likely to exist, if the ground state of matter at large atomic 
number is in the form of a quark fluid, which would then be composed of about equal 
numbers of up, down and strange quarks together with electrons, which give overall 
charge neutrality [37, 95]. The strangeness per unit baryon number is ~ —1. The 
first relativistic models of stars composed of quark matter were computed by Ipser, 
Kislinger and Morley [155] and by Brecher and Caporaso [49], while the first extensive 
study of strange quark star properties is due to Witten [327]. 

The strange quark matter equation of state can be represented by the following 
Unear relation between pressure and energy density 

P = a(e-eo), (27) 

where eo is the energy density at the surface of a bare strange star (neglecting a possible 
thin crust of normal matter). The MIT bag model of strange quark matter involves three 
parameters, the bag constant, ® = eo/4, the mass of the strange quark, ot.s, and the QCD 
coupUng constant, 0.c. The constant a in (27) is equal to 1/3 if one neglects the mass 
of the strange quark, while it takes the value of a = 0.289 for = 250 MeV. When 
measured in units of «60 = '2 / (60 MeV fm~^), the constant B is restricted to be in the 
range 

0.982K Sgo < 1-525, (28) 

assuming rris = 0. The lower limit is set by the requirement of stability of neutrons with 
respect to a spontaneous fusion into strangelets, while the upper limit is determined by 
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the energy per baryon of ^^Fe at zero pressure (930.4 MeV). For other values of nis the 
above limits are modified somewhat. 

A more recent attempt at describing deconfined strange quark matter is the Dey et 
al. EOS [84], which has asymptotic freedom built in. It describes deconfined quarks at 
high densities and confinement at zero pressure. The Dey et al. EOS can be approx- 
imated by a linear relation of the same form as the MIT bag-model strange star EOS 
(27). In such a linear approximation, typical values of the constant a are 0.45 — 0.46 
[124]. 

• Going further. A review of strange quark star properties can be found in [322]. 
Hybrid stars, that have a mixed-phase region of quark and hadronic matter, have 
also been proposed, see e.g. [122]. A study of the relaxation effect in dissipative 
relativistic fluid theories is presented in [197]. 

2.7 Numerical Schemes 

AH available methods for solving the system of equations describing the equilibrium 
of rotating relativistic stars are numerical, as no analytical self-consistent solution for 
both the interior and exterior spacetime has been found. The first numerical solutions 
were obtained by Wilson [325] and by Bonazzola & Schneider [47]. Here, we will 
review the following methods: Hartle's slow rotation formalism, the Newton-Raphson 
linearization scheme due to Butterworth & Ipser [54], a scheme using Green's functions 
by Komatsu et al. [181, 182], a minimal surface scheme due to Neugebauer & Herold 
[234], and spectral-method schemes by Bonazzola et al. [46, 45] and by Ansorg et 
al. [17]. Below we give a description of each method and its various implementations 
(codes). 

2.7.1 Hartle 

To order 0{Q?) the structure of a star changes only by quadrupole terms and the equi- 
librium equations become a set of ordinary differential equations. Hartle's [140, 145] 
method computes rotating stars in this slow-rotation approximation and a review of 
slowly rotating models has been compiled by Datta [80]. Weber et al. [319], [321] also 
implement Hartle's formalism to explore the rotational properties of four new EOSs. 

Weber and Glendenning [320] improve on Hartle's formalism in order to obtain 
a more accurate estimate of the angular velocity at the mass-shedding limit, but their 
models still show large discrepancies compared to corresponding models computed 
without the slow-rotation approximation [262]. Thus, Hartle's formalism is appropriate 
for typical pulsar (and most millisecond pulsar) rotational periods, but it is not the 
method of choice for computing models of rapidly rotating relativistic stars near the 
mass-shedding limit. 
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2.7.2 Butterworth and Ipser (BI) 

The BI scheme [54] solves the four field equations following a Newton-Raphson-like 
linearization and iteration procedure. One starts with a nonrotating model and increases 
the angular velocity in small steps, treating a new rotating model as a linear pertur- 
bation of the previously computed rotating model. Each hnearized field equation is 
discretized and the resulting linear system is solved. The four field equations and the 
hydrostationary equiUbrium equation are solved separately and iteratively until conver- 
gence is achieved. 

Space is truncated at a finite distance from the star and the boundary conditions 
there are imposed by expanding the metric potentials in powers of 1 jr. Angular deriva- 
tives are approximated by high-accuracy formulae and models with density disconti- 
nuities are treated specially at the surface. An equiUbrium model is specified by fixing 
its rest mass and angular velocity. 

The original BI code was used to construct uniform density models and polytropic 
models [54, 53]. Friedman et al. [110, 111], extend the BI code to obtain a large 
number of rapidly rotating models based on a variety of realistic EOSs. Lattimer et al. 
[193] used a code which was also based on the BI scheme to construct rotating stars 
using "exotic" and schematic EOSs, including pion or Kaon condensation and strange 
quark matter. 

2.7.3 Komatsu, Eriguchi and Hachisu (KEH) 

In the KEH scheme [181, 182], the same set of field equations as in BI is used, but the 
three elliptic-type field equations are converted into integral equations using appropri- 
ate Green's functions. The boundary conditions at large distance from the star are thus 
incorporated into the integral equations, but the region of integration is truncated at a 
finite distance from the star. The fourth field equation is an ordinary first-order differ- 
ential equation. The field equations and the equation of hydrostationary equilibrium 
are solved iteratively, fixing the maximum energy density and the ratio of the polar 
radius to the equatorial radius, until convergence is achieved. In [181, 182, 92] the 
original KEH code is used to construct uniformly and differentially rotating stars for 
both polytropic and realistic EOSs. 

Cook, Shapiro and Teukolsky (CST) improve on the KEH scheme by introducing 
a new radial variable which maps the semi-infinite region [0, °°) to the closed region 
[0, 1]. In this way, the region of integration is not truncated and the model converges to 
a higher accuracy. Details of the code are presented in [67] and polytropic and reaUstic 
models are computed in [68] and [69]. 

Stergioulas and Friedman (SF) implement their own KEH code following the CST 
scheme. They improve on the accuracy of the code by a special treatment of the second 
order radial derivative that appears in the source term of the first-order differential 
equation for one of the metric functions. This derivative was introducing a numerical 
error of 1% — 2% in the bulk properties of the most rapidly rotating stars computed in 
the original implementation of the KEH scheme. The SF code is presented in [296] and 
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in [294]. It is available as a public domain code, named rns, and can be downloaded 
from [293]. 

2.7.4 Bonazzola et al. (BGSM) 

In the BGSM scheme [46], the field equations are derived in the 3 + 1 formulation. 
All four chosen equations that describe the gravitational field are of elliptic type. This 
avoids the problem with the second-order radial derivative in the source term of the 
ODE used in BI and KEH. The equations are solved using a spectral method, i.e. all 
functions are expanded in terms of trigonometric functions in both the angular and 
radial directions and a Fast Fourier Transform (FFT) is used to obtain coefficients. 
Outside the star a redefined radial variable is used, which maps infinity to a finite 
distance. 

In [262, 263] the code is used to construct a large number of models based on recent 
EOSs. The accuracy of the computed models is estimated using two general relativistic 
Virial identities, valid for general asymptotically flat spacetimes [128, 42] (see Section 
2.7.7). 

While the field equations used in the BI and KEH schemes assume a perfect fluid, 
isotropic stress-energy tensor, the BGSM formulation makes no assumption about the 
isotropy of Tab- Thus, the BGSM code can compute stars with a magnetic field, a solid 
crust or solid interior and it can also be used to construct rotating boson stars. 

2.7.5 Lorene/rotstar 

The BGSM spectral method has been improved by Bonazzola et al. [45] allowing for 

several domains of integration. One of the domain boundaries is chosen to coincide 
with the surface of the star and a regularization procedure is introduced for the di- 
vergent derivatives at the surface (that appear in the density field when stiff equations 
of state are used). This allows models to be computed that are nearly free of Gibbs 
phenomena at the surface. The same method is also suitable for constructing quasi- 
stationary models of binary neutron stars. The new method has been used in [129] 
for computing models of rapidly rotating strange stars and it has also been used in 
3D computations of the onset of the viscosity-driven instability to bar-mode formation 
[125]. 

2.7.6 Ansorg et al. (AKM) 

A new multi-domain spectral method has been introduced in [17, 16]. The method can 
use several domains inside the star, one for each possible phase transition. Surface- 
adapted coordinates are used and approximated by a two-dimensional Chebyshev- 
expansion. Requiring transition conditions to be satisfied at the boundary of each do- 
main, the field and fluid equations are solved as a free boundary value problem by a 
Newton-Raphson method, starting from an initial guess. The field equations are sim- 
plified by using a corotating reference frame. Applying this new method to the com- 
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putation of rapidly rotating homogeneous relativistic stars, Ansorg et al. achieve near 
machine accuracy, except for configurations at the mass-shedding Umit (see Section 
2.7.8)! The code has been used in a systematic study of uniformly rotating homoge- 
neous stars in general relativity [266]. 

2.7.7 The Virial Identities 

EquiUbrium configurations in Newtonian gravity satisfy the well-known virial relation 

2T + 3{r-l)U + W = 0. (29) 

This can be used to check the accuracy of computed numerical solutions. In general 
relativity, a different identity, valid for a stationary and axisymmetric spacetime was 
found in [39]. More recently, two relativistic virial identities, valid for general asymp- 
totically flat spacetimes, have been derived by Bonazzola and Gourgoulhon [128, 42]. 
The 3-dimensional virial identity (GRV3) [128] is the extension of the Newtonian virial 
identity (29) to general relativity. The 2-dimensional (GRV2) [42] virial identity is the 
generalization of the identity found in [39] (for axisymmetric spacetimes) to general 
asymptotically flat spacetimes. In [42], the Newtonian limit of GRV2, in axisymme- 
try, is also derived. Previously, such a Newtonian identity had only been known for 
spherical configurations [57]. 

The two virial identities are an important tool for checking the accuracy of numer- 
ical models and have been repeatedly used by several authors [46, 262, 263, 235, 17]. 

2.7.8 Direct Comparison of Numerical Codes 

The accuracy of the above numerical codes can be estimated, if one constructs exactly 
the same models with different codes and compares them directly. The first such com- 
parison of rapidly rotating models constructed with the FIP and SF codes is presented 
by Stergioulas and Friedman in [296]. Rapidly rotating models constructed with sev- 
eral EOS's agree to 0.1% — 1.2% in the masses and radii and to better than 2% in any 
other quantity that was compared (angular velocity and momentum, central values of 
metric functions etc.). This is a very satisfactory agreement, considering that the BI 
code was using relatively few grid points, due to limitations of computing power at the 
time of its implementation. 

In [296], it is also shown that a large discrepancy between certain rapidly rotating 
models, constructed with the FIP and KEH codes, that was reported by Eriguchi et al. 
[92], was only due to the fact that a different version of a tabulated EOS was used in 
[92] than by FTP. 

Nozawa et al. [235] have completed an extensive direct comparison of the BGSM, 
SF and the original KEH codes, using a large number of models and equations of state. 
More than twenty different quantities for each model are compared and the relative 
differences range from 10^"^ to 10^^ or better, for smooth equations of state. The 
agreement is excellent for soft polytropes, which shows that all three codes are correct 
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AKM Lorene/ SF SF BGSM KEH 

rotstar (260x400) (70x200) 



Pc 1 



rp/re 


0.7 








le-3 




Q. 


1.41170848318 


9e-6 


3e-4 


3e-3 


le-2 


le-2 


M 


0.135798178809 


2e-4 


2e-5 


2e-3 


9e-3 


2e-2 


Mo 


0.186338658186 


2e-4 


2e-4 


3e-3 


le-2 


2e-3 




0.345476187602 


5e-5 


3e-5 


5e-4 


3e-3 


le-3 


J 


0.0140585992949 


2e-5 


4e-4 


5e-4 


2e-2 


2e-2 




1.70735395213 


le-5 


4e-5 


le-4 


2e-2 


6e-2 




-0.162534082217 


2e-4 


2e-3 


2e-2 


4e-2 


2e-2 


Z^ 


11.3539142587 


7e-6 


7e-5 


le-3 


8e-2 


2e-l 


|GRV3| 


4e-13 


3e-6 


3e-5 


le-3 


4e-3 


le-1 



Table 2: Detailed comparison of the accuracy of different numerical codes in com- 
puting a rapidly rotating, uniform density model. The absolute value of the relative 
error in each quantity, compared to the AKM code, is shown for the numerical codes 
Lorene/rotstar, SF (at two resolutions), BGSM and KEH (see text). The resolutions 
for the SF code are (angular x radial) grid points. See [235] for the definition of the 
various equilibrium quantities. 

and compute the desired models to an accuracy that depends on the number of grid- 
points used to represent the spacetime. 

If one makes the extreme assumption of uniform density, the agreement is at the 
level of 10^^. In the BGSM code this is due to the fact that the spectral expansion 
in terms of trigonometric functions cannot accurately represent functions with discon- 
tinuous first-order derivatives at the surface of the star. In the KEH and SF codes, the 
three-point finite-difference formulae cannot accurately represent derivatives across the 
discontinuous surface of the star. 

The accuracy of the three codes is also estimated by the use of the two Virial iden- 
tities. Overall, the BGSM and SF codes show a better and more consistent agreement 
than the KEH code with BGSM or SF. This is largely due to the fact that the KEH code 
does not integrate over the whole spacetime but within a finite region around the star, 
which introduces some error in the computed models. 

A new direct comparison of different codes is presented by Ansorg et al. [17]. 
Their multi-domain spectral code is compared to the BGSM, KEH and SF codes for a 
particular uniform-density model of a rapidly rotating relativistic star. An extension of 
the detailed comparison in [17], which includes results obtained by the Lorene/rotstar 
code in [125] and by the SF code with higher resolution than the resolution used in 
[235], is shown in Table 2. The comparison confirms that the virial identity GRV3 is a 
good indicator for the accuracy of each code. For the particular model in Table 2, the 



15 



AKM code achieves near double precision machine accuracy, while the Lorene/rotstar 
code has typical relative accuracy of 2 x 10""* to 7 x 10~^ in various quantities. The 
SF code at high resolution comes close to the accuracy of the Lorene/rotstar code for 
this model. Lower accuracies are obtained with the SF, BGSM and KEH codes at the 
resolutions used in [235]. 

The AKM code converges to machine accuracy when a large number of about 24 
expansion coefficients are used at a high computational cost. With significantly fewer 
expansion coefficients (and comparable computational cost to the SF code at high res- 
olution) the achieved accuracy is comparable to the accuracy of the Lorene/rotstar and 
SF codes). Moreover, at the mass-shedding limit, the accuracy of the AKM code re- 
duces to about 5 digits (which is still highly accurate, of course), even with 24 ex- 
pansion coefficients, due to the nonanalytic behaviour of the solution at the surface. 
Nevertheless, the AKM method represents a great achievement, as it is the first method 
to converge to machine accuracy when computing rapidly rotating stars in general rel- 
ativity. 

• Going further. A review of spectral methods in general relativity can be found 
in [4 1] . A formulation for nonaxisymmetric, uniformly rotating equilibrium con- 
figurations in the second post-Newtonian approximation is presented in [20J. 

2.8 Analytic Approximations to the Exterior Spacetime 

The exterior metric of a rapidly rotating neutron star differs considerably from the 
Kerr metric. Only to lowest order in rotation, the two metrics agree [146]. At higher 
order, the multipole moments of the gravitational field created by a rapidly rotating 
compact star are different from the multipole moments of the Kerr field. There have 
been many attempts in the past to find analytic solutions to the Einstein equations 
in the stationary, axisynmietric case, that could describe a rapidly rotating neutron 
star. An interesting solution has been found recently by Manko et al. [218, 219]. 
For non-magnetized sources of zero net charge, the solution reduces to a 3-parameter 
solution, involving the mass, specific angular momentum and a parameter that depends 
on the quadrupole moment of the source. Although it depends explicitly only on the 
quadrupole moment, it approximates the gravitational field of a rapidly rotating star 
with higher non-zero multipole moments. It would be interesting to determine whether 
this analytic quadrupole solution approximates the exterior field of a rapidly rotating 
star more accurately than the quadrupole, 0{Q?'), slow-rotation approximation. 

The above analytic solution (and an earlier one that was not represented in terms 
of rational functions [217]) have been used in studies of energy release during disk 
accretion onto a rapidly rotating neutron star [280, 281]. In [277], a different approxi- 
mation to the exterior spacetime, in the form of a multipole expansion far from the star, 
has been used to derive approximate analytic expressions for the location of the inner- 
most stable circular orbit (ISCO). Even though the analytic solutions in [277] converge 
slowly to an exact numerical solution at the surface of the star, the analytic expressions 
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for the location and angular velocity at the ISCO are in good agreement with numerical 
results. 

2.9 Properties of Equilibrium Models 
2.9.1 Bulk Properties of Equilibrium Models 

Neutron star models constructed with various realistic EOSs have considerably dif- 
ferent bulk properties, due to the large uncertainties in the equation of state at high 
densities. Very compressible (soft) EOSs produce models with small maximum mass, 
small radius, and large rotation rate. On the other hand, less compressible (stiff) EOSs 
produce models with a large maximum mass, large radius, and low rotation rate. 

The gravitational mass, equatorial radius and rotational period of the maximum 
mass model constructed with one of the softest EOSs (EOS B) (1 .63M0, 9.3km, 0.4ms) 
are a factor of two smaller than the mass, radius and period of the corresponding model 
constructed by one of the stiffest EOSs (EOS L) {3.21Mq, 18.3km, 0.8ms). The two 
models differ by a factor of 5 in central energy density and a factor of 8 in the moment 
of inertia! 

Not all properties of the maximum mass models between proposed EOSs differ 

considerably, at least not within groups of similar EOSs. For example, most realistic 
hadronic EOSs predict a maximum mass model with a ratio of rotational to gravita- 
tional energy T /W of 0.11 ±0.02, a dimensionless angular momentum cJ/GM^ of 
0.64 ± 0.06 and an eccentricity of 0.66 ± 0.04, [109]. Hence, between the set of realis- 
tic hadronic EOSs, some properties are directly related to the stiffness of the EOS while 
other properties are rather insensitive to stiffness. On the other hand, if one considers 
strange quark EOSs, then for the maximum mass model, T /W can become a factor of 
about two larger than for hadronic EOSs. 

Compared to nonrotating stars, the effect of rotation is to increase the equatorial 
radius of the star and also to increase the mass that can be sustained at a given central 
energy density. As a result, the mass of the maximum mass rotating model is roughly 
15% — 20% higher than the mass of the maximum mass nonrotating model, for typical 
realistic hadronic EOSs. The corresponding increase in radius is 30% —40%. The 
effect of rotation in increasing the mass and radius becomes more pronounced in the 
case of strange quark EOSs (see Section 2.9.8). 

The deformed shape of a rapidly rotating star creates a distortion, away from spher- 
ical symmetry, in its gravitational field. Far from the star, the dominant multipole mo- 
ment of the rotational distortion is measured by the quadrupole-moment tensor Qab- 
For uniformly rotating, axisymmetric and equatorially symmetric configurations, one 
can define a scalar quadrupole moment Q, which can be extracted from the asymptotic 
expansion, at large r, of the metric function v, as in Equation (10). 

Laarakkers & Poisson [185] numerically compute the scalar quadrupole moment 
Q for several equations of state, using the rotating neutron star code rns [293]. They 
find that for fixed gravitational mass M, the quadrupole moment is given as a simple 
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Figure 1: 2-D surface of equilibrium models for EOS L. The surface is bounded by 
the nonrotating (J — 0) and mass-shedding (Q.~ £Lk) limits and formed by constant J 
and constant Mq sequences (solid lines). The projection of these sequences in the J- 
M plane are shown as long-dashed lines. Also shown are the axisymmetric instability 
sequence (short-dashed line). The projection of the 2-D surface in the J-M plane shows 
an overlapping (see dotted lines). (Figure 7 of Stergioulas and Friedman, ApJ [296].) 



quadratic fit 

e = -^, (30) 

where J is the angular momentum of the star and a is a dimensionless quantity that de- 
pends on the equation of state. The above quadratic fit reproduces Q with a remarkable 
accuracy. The quantity a varies between a 2 for very soft EOSs and a 8 for very 
stiff EOSs, for M = 1.4M© neutron stars. This is considerably different from a Kerr 
black hole, for which a = I [308]. 

For a given zero-temperature EOS, the uniformly rotating equilibrium models form 
a 2-dimensional surface in the 3-dimensional space of central energy density, gravita- 
tional mass and angular momentum [296], as shown in Figure 1 for EOS L. The surface 
is limited by the nonrotating models (/ = 0) and by the models rotating at the mass- 
shedding (Kepler) limit, i.e. at the maximum allowed angular velocity so that the star 
does not shed mass at the equator. Cook et al. [67, 68, 69] have shown that the model 
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with maximum angular velocity does not coincide with the maximum mass model, but 
is generally very close to it in central density and mass. Stergioulas and Friedman [296] 
show that the maximum angular velocity and maximum baryon mass equilibrium mod- 
els are also distinct. The distinction becomes significant in the case where the EOS has 
a large phase transition near the central density of the maximum mass model; otherwise 
the models of maximum mass, baryon mass, angular velocity and angular momentum 
can be considered to coincide for most purposes. 

• Going further. Although rotating relativistic stars are nearly perfectly axisym- 
metric, a small degree of asymmetry (e.g. frozen into the solid crust during its 
formation) can become a source of gravitational waves. A recent review of this 
can be found in [162]. 



2.9.2 Mass-shedding Limit and the Empirical Formula 

Mass-shedding occurs when the angular velocity of the star reaches the angular velocity 
of a particle in a circular Keplerian orbit at the equator, i.e. 

a = ^K, (31) 

where 



c h e 



-I 1/2 

+ C0. (32) 



In differentially rotating stars, even a small amount of differential rotation can increase 
the angular velocity required for mass-shedding significantly. Thus, a newly-born, hot, 
differentially rotating neutron star or a massive compact object created in a binary 
neutron star merger, could be sustained (temporarily) in equilibrium by differential 
rotation, even if a uniformly rotating configuration with same rest mass does not exist. 

In the Newtonian limit the maximum angular velocity of uniformly rotating poly- 
tropic stars is approximately Cl^mx — (2/3)^/^(GM//?^)'/^ (this is derived using the 
Roche model, see [269]). For relativistic stars, the empirical formula [139, 111, 105] 



il^a. =0.67 J (33) 

1/ (U/ 



gives the maximum angular velocity in terms of the mass and radius of the maximum 
mass nonrotating model with an accuracy of 5% — 7%, without actually having to con- 
struct rotating models. A revised empirical formula, using a large set of EOSs, has 
been computed in [138]. 

The empirical formula results from universal proportionality relations that exist 
between the mass and radius of the maximum mass rotating model and those of the 
maximum mass nonrotating model for the same EOS. Lasota et al. [190] find that, for 
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most EOSs, the coefficient in the empirical formula is an almost linear function of the 
parameter 

= (34) 



The Lasota et al. empirical formula 



^max = C{Xs)i , (35) 

y max 

with c{Xs)= 0.468 + 0.378Xs, reproduces the exact values with a relative error of only 
1.5%. 

Weber and Glendenning [319, 320], derive analytically a similar empirical formula 
in the slow rotation approximation. However, the formula they obtain involves the mass 
and radius of the maximum mass rotating configuration, which is different from what 
is involved in (33). 



2.9.3 Upper Limits on Mass and Rotation: Theory vs. Observation 

The maximum mass and minimum period of rotating relativistic stars computed with 
realistic hadronic EOSs from the Amett and Bowers collection [18] are about 3.3M0 
(EOS L) and 0.4 ms (EOS B), while 1.4M0 neutron stars, rotating at the Kepler limit, 
have rotational periods between 0.53 ms (EOS B) and 1 .7 ms (EOS M) [69]. The max- 
imum, accurately measured, neutron star mass is currently still 1 .44M0 (see e.g. [316], 
but there are also indications for 2.0M neutron stars [164]. Core-collapse simulations 
have yielded a bi-modal mass distribution of the remnant, with peaks at about 1.3M0 
and 1.7M0 [311] (the second peak depends on the assumption for the high-density 
EOS - if a soft EOS is assumed, then black hole formation of this mass is implied). 
Compact stars of much higher mass, created in a neutron star binary merger, could be 
temporarily supported against collapse by strong differential rotation [30]. 

When magnetic-field effects are ignored, conservation of angular momentum can 
yield very rapidly rotating neutron stars at birth. Recent simulations of the rotational 
core collapse of evolved rotating progenitors [148, 116] have demonstrated that rota- 
tional core collapse can easily result in the creation of neutron stars with rotational 
periods of the order of 1 ms (and similar initial rotation periods have been estimated 
for neutron stars created in the accretion-induced collapse of a white dwarf [210]). The 
existence of a magnetic field may complicate this picture. Spruit & Phinney [289] 
have presented a model in which a strong internal magnetic field couples the angular 
velocity between core and surface during most evolutionary phases. The core rotation 
decouples from the rotation of the surface only after central carbon depletion takes 
place. Neutron stars born in this way would have very small initial rotation rates, even 
smaller than the ones that have been observed in pulsars associated with supernova 
remnants. In this model, an additional mechanism is required to spin-up the neutron 
star to observed periods. On the other hand, Livio & Pringle [211] argue for a much 
weaker rotational coupling between core and surface by a magnetic field, allowing for 
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the production of more rapidly rotating neutron stars than in [289]. A new investiga- 
tion by Heger et al., yielding intermediate initial rotation rates, is presented in [149]. 
Clearly, more detailed computations are needed to resolve this important question. 

The minimum observed pulsar period is still 1.56ms [184], which is close to the 
experimental sensitivity of most pulsar searches. New pulsar surveys, in principle 
sensitive down to a few tenths of a millisecond, have not been able to detect a sub- 
millisecond pulsar [51, 79, 73, 91]. This not too surprising, as there are several ex- 
planations for the absence of sub-millisecond pulsars. In one model, the minimum 
rotational period of pulsars could be set by the occurrence of the r-mode instability 
in accreting neutron stars in LMXB's [10]. Other models are based on the standard 
magnetospheric model for accretion-induced spin-up [324] or on the idea that gravita- 
tional radiation (produced by accretion-induced quadrupole deformations of the deep 
crust) balances the spin-up torque [35,314]. It has also been suggested [52] that the ab- 
sence of sub-millisecond pulsars in all surveys conducted so far is because they could 
be more likely to be found only in close systems (of orbital period Porb ^ 1 hr), for 
which the current pulsar surveys are still lacking the required sensitivity. The absence 
of sub-millisecond pulsars in wide systems is suggested to be due to the tuming-on of 
the accreting neutron stars as pulsars, in which case the pulsar wind is shown to halt 
further spin-up. 

• Going further. A review by J. L. Friedman on the upper limit on rotation of 
relativistic stars can be found in [107]. 

2.9.4 The Upper Limit on Mass and Rotation Set by Causality 

If one is interested in obtaining upper limits on the mass and rotation rate, indepen- 
dently of the proposed EOSs, one has to rely on fundamental physical principles. In- 
stead of using realistic EOSs, one constructs a set of schematic EOSs that satisfy only a 
minimal set of physical constraints, which represent what we know about the equation 
of state of matter with high confidence. One then searches among all these EOSs to 
obtain the one that gives the maximum mass or minimum period. The minimal set of 
constraints that have been used in such searches are that 

1. the high density EOS matches to the known low density EOS at some matching 

energy density E„, 

2. the matter at high densities satisfies the causaUty constraint (the speed of sound 
is less than the speed of light). 

In relativistic perfect fluids, the speed of sound is the characteristic velocity of the evo- 
lution equations for the fluid and the causality constraint translates into the requirement 

dp/de<l. (36) 

(see [251]). It is assumed that the fluid will stiU behave as a perfect fluid when it is 
perturbed from equilibrium. 
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For nonrotating stars, Rhoades and Ruffini showed that the EOS that satisfies the 
above two constraints and yields the maximum mass consists of a high density region 
as stiff as possible (i.e. at the causal limit, dp/dz = 1), that matches directly to the 
known low density EOS. For a chosen matching density £„, they computed a maximum 
mass of 3.2M0. However, this is not the theoretically maximum mass of nonrotating 
neutron stars, as is often quoted in the literature. Hartle and Sabbadini [143] point out 
that Mmax is sensitive to the matching energy density and Hartle [141] computes M„ax 
as a function of e„ 



M_ = 4.8''^^^«"^^/-^'/^ 



/ ZxlU--gr/cm y/^^^_ (37) 



In the case of rotating stars, Friedman and Ipser [108] assume that the absolute 
maximum mass is obtained by the same EOS as in the nonrotating case and compute 
Mmax as a function of matching density, assuming the BPS EOS holds at low densities. 
A more recent computation [183] uses the FPS EOS at low densities, arriving at a 
similar result as in [108] 

/2xl0iV/cm^i/2 



M' 



/Z X lU""gr/cm \ 1/^ 
•1 Mq, (38) 



where, 2 x lO'^gr/cm^ is roughly nuclear saturation density for the FPS EOS. 

A first estimate of the absolute minimum period of uniformly rotating, gravita- 
tionally bound stars was computed by Glenderming [121] by constructing nonrotating 
models and using the empirical formula (33) to estimate the minimum period. Ko- 
randa, Stergioulas and Friedman [183] improve on Glendenning's results by construct- 
ing accurate rapidly rotating models and show that Glendenning's results are accurate 
to within the accuracy of the empirical formula. 

Furthermore, they show that the EOS satisfying the minimal set of constraints and 
yielding the minimum period star consists of a high density region at the causal limit 
(CL EOS), f = (e — ec), (where ec is the lowest energy density of this region), which 
is matched to the known low density EOS through an intermediate constant pressure 
region (that would correspond to a first-order phase transition). Thus, the EOS yielding 
absolute minimum period models is as stiff as possible at the central density of the star 
(to sustain a large enough mass) and as soft as possible in the crust, in order to have the 
smallest possible radius (and rotational period). 

The absolute minimum period of uniformly rotating stars is an (almost linear) func- 
tion of the maximum observed mass of nonrotating neutron stars 

= 0.28ms + 0.2(M— /Mq - 1 .44)ms, (39) 

and is rather insensitive to the matching density (the above result was computed for 
a matching number density of 0. Ifm"^). In [183], it is also shown that an absolute limit 
on the minimum period exists even without requiring that the EOS matches to a known 
low density EOS, i.e. if the CL EOS P = (e — Ec) terminates at a surface energy density 
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of ec. This is not so for the causal limit on the maximum mass. Thus, without matching 
to a low-density EOS, the causaUty Umit on Pmin is lowered by only 3%, which shows 
that the currently known part of the nuclear EOS plays a negligible role in determining 
the absolute upper limit on the rotation of uniformly rotating, gravitationally bound 
stars. 

The above results have been confirmed in [136], where it is shown that the CL EOS 
has Xs = 0.7081, independent of ec, and the empirical formula (35) reproduces the 
numerical result (39) to within 2%. 

2.9.5 Supramassive Stars and Spin-Up Prior to Collapse 

Since rotation increases the mass that a neutron star of given central density can sup- 
port, there exist sequences of neutron stars with constant baryon mass that have no 
nonrotating member. Such sequences are called supramassive, as opposed to normal 
sequences that do have a nonrotating member. A nonrotating star can become supra- 
massive by accreting matter and spinning-up to large rotation rates; in another sce- 
nario, neutron stars could be born supramassive after a core collapse. A supramassive 
star evolves along a sequence of constant baryon mass, slowly losing angular momen- 
tum. Eventually, the star reaches a point where it becomes unstable to axisymmetric 
perturbations and collapses to a black hole. 

In a neutron star binary merger, prompt collapse to a black hole can be avoided if 
the equation of state is sufficiently stiff and/or the equiUbrium is supported by strong 
differential rotation. The maximum mass of differentially rotating supramassive neu- 
tron stars can be significantly larger than in the case of uniform rotation. A detailed 
study of this mass-increase has recently appeared in [213]. 

Cook et al. [67, 68, 69] have discovered that a supramassive relativistic star ap- 
proaching the axisymmetric instability, will actually spin up before collapse, even 
though it loses angular momentum. This, potentially observable, effect is independent 
of the equation of state and it is more pronounced for rapidly rotating massive stars. 
Similarly, stars can spin up by loss of angular momentum near the mass-shedding limit, 
if the equation of state is extremely stiff or extremely soft. 

If the equation of state features a phase transition to e.g. quark matter, then the spin- 
up region is very large and most milhsecond pulsars (if supramassive) would need to be 
spinning up [290] - the absence of spin-up in known milhsecond pulsars indicates that 
either large phase transitions do not occur, or that the equation of state is sufficiently 
stiff, so that milhsecond pulsars are not supramassive. 

2.9.6 Rotating IMagnetized Neutron Stars 

The presence of a magnetic field has been ignored in the models of rapidly rotating 
relativistic stars that were considered in the previous sections. The reason is that the 
observed surface dipole magnetic field strength of pulsars ranges between 10*^ G and 
2 X lO'^ G. These values of the magnetic field strength imply a magnetic field energy 
density that is too small compared to the energy density of the fluid, to significantly 
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affect the structure of a neutron star. However, one cannot exclude the existence of 
neutron stars with higher magnetic field strengths or the possibiUty that neutron stars 
are born with much stronger magnetic fields, which then decay to the observed values 
(of course, there are also many arguments against magnetic field decay in neutron stars 
[242]). In addition, even though moderate magnetic field strengths do not alter the bulk 
properties of neutron stars, they may have an effect on the damping or growth rate of 
various perturbations of an equilibrium star, affecting its stability. For these reasons, 
a fully relativistic description of magnetized neutron stars is desirable and, in fact, 
Bocquet et al. [36] achieved the first numerical computation of such configurations. 
Here we give a brief summary of their work: 

A magnetized relativistic star in equiUbrium can be described by the coupled Einstein- 
Maxwell field equations for stationary, axisymmetric rotating objects with internal 
electric currents. The stress-energy tensor includes the electromagnetic energy den- 
sity and is non-isotropic (in contrast to the isotropic perfect fluid stress-energy tensor). 
The equilibrium of the matter is given not only by the balance between the gravita- 
tional force, centrifugal force and the pressure gradient, but the Lorentz force due to 
the electric currents also enters the balance. For simplicity, Bocquet et al. consider 
only poloidal magnetic fields, which preserve the circularity of the spacetime. Also, 
they only consider stationary configurations, which excludes magnetic dipole moments 
non-aligned with the rotation axis, since in that case the star emits electromagnetic and 
gravitational waves. The assumption of stationarity implies that the fluid is necessarily 
rigidly rotating (if the matter has infinite conductivity) [46]. Under these assumptions, 
the electromagnetic field tensor F"^ is derived from a potential four- vector Aa with 
only two non-vanishing components. At and A^f, which are solutions of a scalar Pois- 
son and a vector Poisson equation respectively. Thus, the two equations describing the 
electromagnetic field are of similar type as the four field equations that describe the 
gravitational field. 

For magnetic field strengths larger than about lO'"* G, one observes significant ef- 
fects, such as a flattening of the equilibrium configuration. There exists a maximum 
value of the magnetic field strength, of the order of lO'^ G, for which the magnetic 
field pressure at the center of the star equals the fluid pressure. Above this value no 
stationary configuration can exist. 

A strong magnetic field allows a maximum mass configuration with larger M„iax 
than for the same EOS with no magnetic field and this is in analogy with the increase 
of Mmax induced by rotation. For nonrotating stars, the increase in Mmax, due to a strong 
magnetic field, is 13% — 29%, depending on the EOS. Correspondingly, the maximum 
allowed angular velocity, for a given EOS, also increases in the presence of a strong 
magnetic field. 

Another application of general-relativistic E/M theory in neutron stars is the study 
of the evolution of the magnetic field during pulsar spin-down. A detailed analysis of 
the evolution equations of the E/M field in a slowly rotating magnetized neutron star 
has revealed that effects due to the spacetime curvature and due to the rotational frame- 
dragging are present in the induction equations, when one assumes finite electrical 
conductivity (see [253] and references therein). Numerical solutions of the evolution 
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equations of the E/M have shown, however, that for realistic values of the electrical 
conductivity, the above relativistic effects are small, even in the case of rapid rotation 
[338]. 

• Going further. An 0{Q.) slow-rotation approach for the construction of rotating 
magnetized relativistic stars is presented in [134]. 

2.9.7 Rapidly Rotating Proto-Neutron Stars 

Following the gravitational collapse of a massive stellar core, a proto-neutron star 
(PNS) is born. Initially it has a large radius of about 100 km and a temperature of 
50-100 MeV. The PNS may be born with a large rotational kinetic energy and initially 
it will be differentially rotating. Due to the violent nature of the gravitational collapse, 
the PNS pulsates heavily, emitting significant amounts of gravitational radiation. Af- 
ter a few hundred pulsational periods, bulk viscosity will damp the pulsations signifi- 
cantly. Rapid cooling due to deleptonization transforms the PNS into a hot neutron star 
of r ^ 10 MeV, shortly after its formation. In addition, viscosity or other mechanisms 
(see Sec. 2.1) enforce uniform rotation and the neutron star becomes quasi- stationary. 
Since the details of the PNS evolution determine the properties of the resulting cold 
NSs, proto-neutron stars need to be modeled realistically in order to understand the 
structure of cold neutron stars. 

Hashimoto et al. [147] and Goussard et al. [130] construct fully relativistic models 
of rapidly rotating, hot proto-neutron stars. The authors use finite-temperature EOSs 
[238, 192], to model the interior of PNSs. Important (but largely unknown) parameters, 
that determine the local state of matter, are the lepton fraction 7/ and the temperature 
profile. Hashimoto et al. consider only the limiting case of zero lepton fraction F/ = 
and classical isothermality, while Goussard et al. consider several non-zero values for 
Yi and two different limiting temperature profiles - a constant entropy profile and a 
relativistic isothermal profile. In both [147] and [238], differential rotation is neglected 
to a first approximation. 

The construction of numerical models with the above assumptions shows that, due 
to the high temperature and the presence of trapped neutrinos, PNSs have a significantly 
larger radius than cold NSs. These two effects give the PNS an extended envelope 
which, however, contains only roughly 0.1% of the total mass of the star. This outer 
layer cools more rapidly than the interior and becomes transparent to neutrinos, while 
the core of the star remains hot and neutrino opaque for a longer time. The two regions 
are separated by the "neutrino- sphere". 

Compared to the T = case, an isothermal EOS with temperature of 25MeV has 
a maximum mass model of only slightly larger mass. In contrast, an isentropic EOS 
with a nonzero trapped lepton number features a maximum mass model that has a 
considerably lower mass than the corresponding model in the T = case and a stable 
PNS transforms to a stable neutron star. If, however, one considers the hypothetical 
case of a large amplitude phase transition which softens the cold EOS (such as a Kaon 
condensate), then Mmax of cold neutron stars is lower than M^ax of PNSs and a stable 
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Figure 2: Iso-energy density lines of a differentially rotating proto-neutron star at the 
mass-shedding limit, of rest mass Mq =1.5 Mq. (Figure 5a of Goussard, Haensel and 
Zdunik, A&A [131]; used with permission.) 



PNS with maximum mass will collapse to a black hole after the initial cooling period. 
This scenario of delayed collapse of nascent neutron stars has been proposed by Brown 
and Bethe [50] and investigated by Baumgarte et al. [31]. 

An analysis of radial stability of PNSs [123] shows that, for hot PNSs, the max- 
imum angular velocity model almost coincides with the maximum mass model, as is 
also the case for cold EOSs. 

Because of their increased radius, PNSs have a different mass-shedding limit than 
cold NSs. For an isothermal profile, the mass-shedding hmit proves to be sensitive to 
the exact location of the neutrino sphere. For the EOSs considered in [147] and [130] 
PNSs have a maximum angular velocity that is considerably less than the maximum 
angular velocity allowed by the cold EOS. Stars that have nonrotating counterparts 
(i.e. that belong to a normal sequence) contract and speed up while they cool down. 
The final star with maximum rotation is thus closer to the mass-shedding limit of cold 
stars than was the hot PNS with maximum rotation. Surprisingly, stars belonging to 
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a supramassive sequence exhibit the opposite behavior. If one assumes that a PNS 
evolves without losing angular momentum or accreting mass, then a cold neutron star 
produced by the cooling of a hot PNS has a smaller angular velocity than its progenitor. 
This purely relativistic effect was pointed out in [147] and confirmed in [130]. 

It should be noted here, that a small amount of differential rotation significantly 
affects the mass-shedding limit, allowing more massive stars to exist than uniform 
rotation allows. Taking differential rotation into account, Goussard et al. [131] suggest 
that proto-neutron stars created in a gravitational collapse cannot spin faster than 1.7 
ms. A similar result has been obtained by Strobel et al. [299]. The structure of a 
differentially rotating proto-neutron star at the mass-shedding limit is shown in Figure 
2. The outer layers of the star form an extended disk-like structure. 

The above stringent limits on the initial period of neutron stars are obtained assum- 
ing that the PNS evolves in a quasi- stationary manner along a sequence of equilibrium 
models. It is not clear whether these limits will remain valid, if one studies the early 
evolution of PNS without the above assumption. It is conceivable that the thin hot en- 
velope surrounding the PNS does not affect the dynamics of the bulk of the star. If the 
bulk of the star rotates faster than the (stationary) mass-shedding limit of a PNS model, 
then the hot envelope will simply be shed away from the star in the equatorial region, 
if it cannot remain bounded to the star even when differentially rotating. Such a fully 
dynamical study is needed to obtain an accurate upper limit on the rotation of neutron 
stars. 

• Going further. The thermal history and evolutionary tracks of rotating PNS (in 
the 2nd-order slow-rotation approximation) have been studied recently in [301]. 

2.9.8 Rotating Strange Quark Stars 

Most rotational properties of strange quark stars differ considerably from the proper- 
ties of rotating stars constructed with hadronic EOSs. First models of rapidly rotating 
strange quark stars were computed by Friedman ([106], quoted in [120, 119]) and by 
Lattimer et al. [193]. Colpi & Miller [64] use the 0{Q?') approximation and find that 
the spin of strange stars (newly-born, or spun-up by accretion) may be limited by the 
CFS-instability to the I = m = 2 /-mode, since rapidly rotating strange stars tend to 
have T/W > 0.14. Rapidly rotating strange stars at the mass-shedding limit have been 
computed first by Gourgoulhon et al. [129] and the structure of a representative model 
is displayed in Figure 3. 

Nonrotating strange stars obey scaling relations with the constant ® in the MIT bag- 
model of the strange quark matter EOS (Section 2.6.3) and Gourgoulhon et al. [129] 
also obtain scaling relations for the model with maximum rotation rate. The maximum 
angular velocity scales as 



while the allowed range of S impUes an allowed range of 0.5 13 ms < Pmin < 0.640 ms. 
The empirical formula (33) also holds for rotating strange stars with an accuracy of bet- 
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Figure 3: Meridional plane cross section of a rapidly rotating strange star at the mass- 
shedding limit, obtained with a multi-domain spectral code. The various lines are 
isocontours of the log-enthalpy H, as defined in [129]. Solid lines indicate a positive 
value ofH and dashed lines a negative value ( vacuum). The thick solid line denotes the 
stellar surface. The thick dot-dashed line denotes the boundary between the two com- 
putational domains. (Figure 4 of Gourgoulhon, Haensel, Livine, Paluch, Bonazzola 
and March in A&A [129]; used with permission.) 
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ter than 2%. A derivation of the empirical formula in the case of strange stars, starting 
from first principles, has been presented by Cheng & Harko [60], who found that some 
properties of rapidly rotating strange stars can be reproduced by approximating the 
exterior spacetime by the Kerr metric. 

Since both the maximum mass nonrotating and maximum mass rotating models 
obey similar scalings with B, the ratios 

^ J 44 ^nax ^ J 54 ^4j) 

ii^stat ' Dstat ' ^ ' 

^"max ^^max 

are independent of ® (where /^max is the radius of the maximum mass model). The 
maximum mass increases by 44% and the radius of the maximum mass model by 54%, 
while the corresponding increase for hadronic stars is, at best, ~ 20% and ~ 40%, 
correspondingly. The rotational properties of strange star models which are based on 
the Dey et al. EOS [84] are similar to those of the MIT bag-model EOS [37, 327, 95], 
but some quantitative differences exist [124]. 

Accreting strange stars in LMXBs will follow different evolutionary paths in a 
mass vs. central energy density diagram, than accreting hadronic stars [343]. When 
(and if) strange stars reach the mass-shedding limit, the ISCO still exists [298] (while it 
disappears for most hadronic EOSs). Stergioulas, Kluzniak and Bulik [298] show that 
the radius and location of the ISCO for the sequence of mass-shedding models also 
scales as while the angular velocity of particles in circular orbit at the ISCO 

scales as « ^1^. Additional scalings with the constant a in the strange quark EOS (that 
were proposed in [193]) are found to hold within an accuracy of better than ~ 9% for 
the mass-shedding sequence 

Roca^/"^, iloca-i/8. (42) 

In addition, it is found that models at mass-shedding can have T /W as large as 0.28 for 
M= 1.34M0. 

As strange quark stars are very compact, the angular velocity at the ISCO can 
become very large. If the 1066 Hz upper QPO frequency in 4U 1820-30 (see [165] 
and references therein) is the frequency at the ISCO, then it rules out most models of 
slowly-rotating strange stars in LMXBs. However, in [298] it is shown that rapidly 
rotating bare strange stars are still compatible with this observation, as they can have 
ISCO frequencies < 1 kHz even for 1 .4 M© models. On the other hand, if strange stars 
have a thin solid crust, the ISCO frequency at the mass-shedding limit increases by 
about 10% (compared to a bare strange star of the same mass) and the above observa- 
tional requirement is only satisfied for slowly rotating models near the maximum non- 
rotating mass, assuming some specific values of the parameters in the strange star EOS 
[344, 342]. Moderately rotating strange stars, with spin frequencies around 300Hz 
can also be accommodated for some values of the coupling constant etc [340] (see 
also [127] for a detailed study of the ISCO frequency for rotating strange stars). The 
1066 Hz requirement for the ISCO frequency depends, of course, on the adopted model 
of kHz QPOs in LMXBs and other models exist (see next section). 
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If strange stars can have a solid crust, then the density at the bottom of the crust is 
the neutron drip density End 4.1 x 10^^ g cm~^, as neutrons are absorbed by strange 
quark matter. A strong electric field separates the nuclei of the crust from the quark 
plasma. In general, the mass of the crust that a strange star can support is very small, of 
the order of 10~^ Mq . Rapid rotation increases by a few times the mass of the crust and 
the thickness at the equator becomes much larger than the thickness at the poles [342]. 
Zdunik, Haensel & Gourgoulhon [342] also find that the mass Mcmst and thickness fcrust 
of the crust can be expanded in powers of the spin frequency V3 = v/(10^ Hz) as 



where a subscript "0" denotes nonrotating values. For v < 500 Hz, the above expansion 
agrees well with the quadratic expansion derived previously by Glendenning & Weber 
[118]. In a spinning down magnetized strange quark star with crust, part of the crust 
will gradually dissolve into strange quark matter, in a strongly exothermic process. In 
[342], it is estimated that the heating due to deconfinement may exceed the neutrino 
luminosity from the core of a strange star older than ^ 1000 years and may therefore 
influence the cooling of this compact object (see also [336]). 

2.10 Rotating Relativistic Stars in mXBs 
2.10.1 Particle orbits and kHz quasi-periodic oscillations 

In the last few years. X-ray observations of accreting sources in Low Mass X-ray Bina- 
ries (LMXBs) have revealed a rich phenomenology that is waiting to be interpreted cor- 
rectly and could lead to significant advances in our understanding of compact objects 
(see [189, 315, 249]). The most important feature of these sources is the observation of 
(in most cases) twin kHz quasi-periodic oscillations (QPOs) . The high frequency of 
these variabilities and their quasi-periodic nature are evidence that they are produced 
in high-velocity flows near the surface of the compact star. To date, there exist a large 
number of different theoretical models that attempt to explain the origin of these oscil- 
lations. No consensus has been reached, yet, but once a credible explanation is found, 
it will lead to important constraints on the properties of the compact object that is the 
source of the gravitational field in which the kHz oscillations take place. The compact 
stars in LMXBs are spun up by accretion, so that many of them may be rotating rapidly 
and the correct inclusion of rotational effects in the theoretical models for kHz QPOs is 
important. Under simplifying assumptions for the angular momentum and mass evolu- 
tion during accretion, one can use accurate rapidly rotating relativistic models to follow 
the possible evolutionary tracks of compact stars in LMXBs [70, 343]. 

In most theoretical models, one or both kHz QPO frequencies are associated with 
the orbital motion of inhomogeneities or blobs in a thin accretion disk. In the actual 
calculations, the frequencies are computed in the approximation of an orbiting test 
particle, neglecting pressure terms. For most equations of state, stars that are massive 



Mcrust = Mcrust,0(l+0.24V3-^0.16V3), 

W = W,o(l+0.4vf-F0.3vf), 



(43) 
(44) 
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enough posses an innermost-stable circular orbit (ISCO) and the orbital frequency at 
the ISCO has been proposed to be one of the two observed frequencies. To first order in 
the rotation rate, the orbital frequency at the prograde ISCO is given by (see Kluzniak, 
Michelson and Wagoner [168J) 



where j = J /M^. At larger rotation rates, higher order contributions of j as well as 
contributions from the quadrupole moment Q become important and an approximate 
expression has been derived by Shibata and Sasaki [277] , which, when written as above 
and truncated to the lowest order contribution of Q and to 0{f), becomes 



where Q2 = -QjM^. 

Notice that, while rotation increases the orbital frequency at the ISCO, the quadrupole 
moment has the opposite effect, which can become important for rapidly rotating mod- 
els. Numerical evaluations of /isco for rapidly rotating stars have been used in [227] 
to arrive at constraints on the properties of the accreting compact object. 

In other models, orbits of particles that are eccentric and slightly tilted with respect 
to the equatorial plane are involved. For eccentric orbits, the periastron advances with a 
frequency V/>a that is the difference between the Keplerian frequency of azimuthal mo- 
tion Va: and the radial epicyclic frequency Vr. On the other hand, particles in slightly 
tilted orbits fail to return to the initial displacement \|/ from the equatorial plane, after 
a full revolution around the star. This introduces a nodal precession frequency V/>4, 
which is the difference between V/f and the frequency of the motion out of the orbital 
plane (meridional frequency) Vy. Explicit expressions for the above frequencies, in 
the gravitational field of a rapidly rotating neutron star, have been derived recently by 
Markovic [220], while in [221] highly eccentric orbits are considered. Morsink and 
Stella [229] compute the nodal precession frequency for a wide range of neutron star 
masses and equations of state and (in a post-Newtonian analysis) separate the preces- 
sion caused by the Lense-Thirring (frame-dragging) effect from the precession caused 
by the quadrupole moment of the star. The nodal and periastron precession of inclined 
orbits have also been studied using an approximate analytic solution for the exterior 
gravitational field of rapidly rotating stars [279]. These precession frequencies are 
relativistic effects and have been used in several models to explain the kHz QPO fre- 
quencies [292, 250, 1, 166, 3]. 

It is worth mentioning that it has recently been found that an ISCO also exists in 
Newtonian gravity, for models of rapidly rotating low-mass strange stars. The instabil- 
ity in the circular orbits is produced by the large oblateness of the star [167, 341, 3]. 
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2.10.2 Angular Momentum Conservation During Burst Oscillations 

Some sources in LMXBs show signatures of type I X-ray bursts, which are thermonu- 
clear flashes on the surface of the compact star [195]. Such bursts show nearly-coherent 
oscillations in the range 270-620 Hz (see [315, 300] for recent reviews). One interpre- 
tation of the burst oscillations is that they are the result of rotational modulation of 
surface asymmetries during the burst. In such case, the oscillation frequency should 
be nearly equal to the spin frequency of the star This model currently has difficulties 
in explaining some observed properties, such as the oscillations seen in the tail of the 
burst, the frequency increase during the burst, the need for two anti-podal hot spots in 
some sources that ignite at the same time etc. Alternative models also exist (see e.g. 



In the spin-frequency interpretation, the increase in the oscillation frequency by 
a few Hz during the burst is explained as follows: The burning shell is supposed to 
first decouple from the neutron star and then gradually settle down onto the surface. 
By angular momentum conservation, the shell spins-up, giving rise to the observed 
frequency increase. Gumming et al. [74] compute the expected spin-up in full general 
relativity and taking into account rapid rotation. Assuming that the angular momentum 
per unit mass is conserved, the change in angular velocity with radius is given by 



where R is the equatorial radius of the star and all quantities are evaluated at the equator. 

The slow-rotation limit of the above result was derived previously by Abramowicz et 
al. [1]. The fractional change in angular velocity during spin-up can then be estimated 
as 



where Ar is the coordinate expansion of the burning shell, a quantity that depends 
on the shell's composition. Gumming et al. find that the spin down expected if the 
atmosphere rotates rigidly is a factor of two to three times smaller than observed values. 
More detailed modeling is needed to fully explain the origin and properties of burst 
osciUations. 

• Going further. A very interesting topic is the modeling of the expected X-ray 

spectrum of an accretion disk in the gravitational field of a rapidly rotating neu- 
tron star as it could lead to observational constraints on the source of the gravi- 
tational field, see e.g. [304, 280, 281, 34, 33], where work initiated by Kluzniak 
and Wilson [169] in the slow-rotation limit is extended to rapidly rotating rela- 
tivistic stars. 



[248]). 
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3 Oscillations and Stability 



The study of oscillations of relativistic stars is motivated by the prospect of detecting 
such oscillations in electromagnetic or gravitational wave signals. In the same way that 
heUoseismology is providing us with information about the interior of the Sun, the ob- 
servational identification of oscillation frequencies of relativistic stars could constrain 
the high-density equation of state [11]. The oscillations could be excited after a core 
collapse or during the final stages of a neutron star binary merger. Rapidly rotating 
relativistic stars can become unstable to the emission of gravitational waves. 

When the oscillations of an equilibrium star are of small magnitude compared to its 
radius, it will suffice to approximate them as linear perturbations. Such perturbations 
can be described in two equivalent ways. In the Lagrangian approach, one studies the 
changes in a given fluid element as it oscillates about its equilibrium position. In the 
Eulerian approach, one studies the change in fluid variables at a fixed point in space. 
Both approaches have their strengths and weaknesses. 

In the Newtonian limit, the Lagrangian approach has been used to develop vari- 
ational principles [214, 115] but the Eulerian approach proved to be more suitable 
for numerical computations of mode frequencies and eigenfunctions [159, 216, 156, 
158, 157]. Clement [62] used the Lagrangian approach to obtain axisymmetric normal 
modes of rotating stars, while nonaxisymmetric solutions were obtained in the La- 
grangian approach by Imamura et al. [153] and in the Eulerian approach by Managan 
[216] and Ipser and Lindblom [156]. While a lot has been learned from Newtonian 
studies, in the following we wiU focus on the relativistic treatment of osciUations of 
rotating stars. 

3.1 Quasi-Normal IVIodes of Oscillation 

A general linear perturbation of the energy density in a static and spherically symmetric 
relativistic star can be written as a sum of quasi-normal modes that are characterized 
by the indices {l,m) of the spherical harmonic functions and have angular and time- 
dependence of the form 

5£ - f{r)Pp (cos e)e' W+M ^ (49) 

where 6 indicates the Eulerian perturbation of a quantity, co, is the angular frequency 
of the mode, as measured by a distant inertial observer, /(r) represents the radial de- 
pendence of the perturbation and P^{cosQ) are the associated Legendre polynomials. 
Normal modes of nonrotating stars are degenerate in m and it suffices to study the 
axisymmetric (m = 0) case. 

The Eulerian perturbation in the fluid 4- velocity, 5m", can be expressed in terms 
of vector harmonics, while the metric perturbation, 5gab, can be expressed in terms of 
spherical, vector and tensor harmonics. These are either of "polar" or "axial" parity. 
Here, parity is defined to be the change in sign under a combination of reflection in 
the equatorial plane and rotation by 7t. A polar perturbation has parity (— 1)^ while an 
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axial perturbation has parity (—1)'+'. Because of the spherical background, the polar 
and axial perturbations of a nonrotating star are completely decoupled. 

A normal mode solution satisfies the perturbed gravitational field equations 

5{G"'' - STtT"^) = 0, (50) 

and the perturbation of the conservation of the stress-energy tensor 

5(V„r^)=0, (51) 

with suitable boundary conditions at the center of the star and at infinity. The latter 
equation is decomposed into an equation for the perturbation in the energy density 
5e and into equations for the three spatial components of the perturbation in the 4- 
velocity, 6m". As Unear perturbation have a gauge freedom, at most six components of 
the perturbed field equations (50) need to be considered. 

For a given pair {l,m),a solution exists for any value of the frequency co„ consisting 
of a mixture of ingoing- and outgoing-wave parts. Outgoing quasi-normal modes are 
defined by the discrete set of eigenfrequencies for which there are no incoming waves 
at infinity. These are the modes that will be excited in various astrophysical situations. 

The main modes of pulsation that are known to exist in relativistic stars have been 
classified as follows (/o and Xo are typical frequencies and damping times of the most 
important modes in the nonrotating Umit): 

1. Polar fluid modes 

Are slowly damped modes analogous to the Newtonian fluid pulsations: 

• /(undamental)-mode: surface mode due to the interface between the star 

and its surroundings (/o ~ 2kHz, Xq < Is), 

• jp(ressure)-modes: nearly radial (/o > 4kHz, Tq > Is), 

• ^(ravity modes): nearly tangential, only exist in stars that are non-isentropic 
or that have a composition gradient or first-order phase transition (/o < 
500Hz, To > 5s). 

2. Axial and Hybrid fluid modes 

• inertial modes: degenerate at zero frequency in nonrotating stars. In a rotat- 
ing star, some inertial modes are generically unstable to the CFS-instability; 
frequencies from zero to kHz, growth times inversely proportional to a high 
power of the star's angular velocity. Hybrid inertial modes have both axial 
and polar parts even in the limit of no rotation. 

• r(otation) modes: a special case of inertial modes that reduce to the clas- 
sical axial r-modes in the Newtonian limit. Generically unstable to the 
CFS-instabiUty with growth times as short as a few seconds at high rota- 
tion rates. 
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3. Polar and axial spacetime modes 

• w(ave)-modes: Analogous to the quasi-normal modes of a black hole (very 
weak coupling to the fluid). High frequency, strongly damped modes (/o > 
6kHz, To ~ 0.1 ms). 

For a more detailed description of various types of osciUation modes see [176, 175, 
223, 55, 174]. 

3.2 Effect of Rotation on Quasi-Normal Modes 

In a continuous sequence of rotating stars that includes a nonrotating member, a quasi- 
normal mode of index / is defined as the mode which, in the nonrotating limit, reduces 
to the quasi-normal mode of the same index /. Rotation has several effects on the modes 
of a corresponding nonrotating star: 

1. The degeneracy in the index m is removed and a nonrotating mode of index / is 
split into 2/ -|- 1 different {l,m) modes. 

2. Prograde (m < 0) modes are now different from retrograde {m > 0) modes. 

3. A rotating "polar" Z-mode consists of a sum of purely polar and purely axial 
terms [294], e.g. for / = m 

oo 

L(P/+2;'+Az+2;'±i), (52) 

l'=0 

that is, rotation couples a polar /-term to an axial Z ± 1 term (the coupling to the 
l + l term is, however, strongly favoured over the couphng to the / — 1 term [59]). 
Similarly, for a rotating "axial" mode with l=m 

oo 

AT'-^I,iAi+2V+Pi+2V±i), (53) 
i'=o 

4. Frequencies and damping times are shifted. In general, frequencies (in the iner- 
tial frame) of prograde modes increase, while those of retrograde modes decrease 
with increasing rate of rotation. 

5. In rapidly rotating stars, apparent intersections between higher-order modes of 
different I can occur. In such cases, the shape of the eigenfunction is used in the 
mode classification. 

In rotating stars, quasi-normal modes of oscillation have only been studied in the 
slow-rotation limit, in the post-Newtonian and in the Cowling Approximations. The 
solution of the fully-relativistic perturbation equations for a rapidly rotating star is stiU 
a very challenging task and only recently they have been solved for zero- frequency 
(neutral) modes [294, 297]. First frequencies of quasi-radial modes have now been 
obtained through 3D numerical time-evolutions of the nonlinear equations [102]. 
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• Going further. The equations that describe oscillations of the solid crust of a 
rapidly rotating relativistic star are derived by Priou in [246]. The effects of 
superfluid hydrodynamics on the oscillations of neutron stars have been investi- 
gated by several authors, see e.g. [202, 65, 6, 8] and references therein. 

3.3 Axisymmetric Perturbations 

3.3.1 Secular and Dynamical Axisynunetric Instability 

Along a sequence of nonrotating relativistic stars with increasing central energy den- 
sity, there is always a model for which the mass becomes maximum. The maximum- 
mass turning point marks the onset of an instability in the fundamental radial pulsation 
mode of the star. 

Applying the turning point theorem provided by Sorkin [287], Friedman, Ipser and 
Sorkin [1 12] show that in the case of rotating stars a secular axisymmetric instability 
sets in when the mass becomes maximum along a sequence of constant angular mo- 
mentum. An equivalent criterion (implied in [112]) is provided by Cook et al. [67] 
: the secular axisymmetric instability sets in when the angular momentum becomes 
minimum along a sequence of constant rest mass. The instability first develops on a 
secular timescale that is set by the time required for viscosity to redistribute the star's 
angular momentum. This timescale is long compared to the dynamical timescale and 
comparable to the spin-up time following a pulsar glitch. Eventually, the star encoun- 
ters the onset of dynamical instability and collapses to a black hole (see [275] for recent 
numerical simulations). Thus, the onset of the secular instability to axisymmetric per- 
turbations separates stable neutron stars from neutron stars that will collapse to a black 
hole. 

Goussard et al. [130] extend the stability criterion to hot protoneutron stars with 
nonzero total entropy. In this case, the loss of stabiUty is marked by the configuration 
with minimum angular momentum along a sequence of both constant rest mass and 
total entropy. In the nonrotating limit, Gondek et al. [123] compute frequencies and 
eigenfunctions of radial pulsations of hot proto-neutron stars and verify that the secular 
instability sets in at the maximum mass turning point, as is the case for cold neutron 
stars. 

3.3.2 Axisymmetric Pulsation Modes 

Axisymmetric (m = 0) pulsations in rotating relativistic stars could be excited in a num- 
ber of different astrophysical scenarios, such as during core collapse, in star quakes 
induced by the secular spin-down of a pulsar or during a large phase transition, in the 
merger of two relativistic stars in a binary system etc. Due to rotational couplings, the 
eigenfunction of any axisymmetric mode will involve a sum of various spherical har- 
monics F^", so that even the quasi-radial modes (with lowest-order 1 = contribution) 
would, in principle, radiate gravitational waves. 
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Figure 4: Apparent intersection (due to avoided crossing) of the axisymmetric first 

quasi-radial overtone (Hi) and the first overtone of the 1=4, p-mode (in the Cowl- 
ing approximation). Frequencies are normalized by Pc/47t, where Pc is the central 
energy density of the star. The rotational frequency frot at the mass-shedding limit is 
0.597 (in the above units). Along continuous sequences of computed frequencies, mode 
eigenfunctions are exchanged at the avoided crossing. Defining quasi-normal mode se- 
quences by the shape of their eigenfunction, the Hi sequence (filled boxes) appears to 
intersect with the '^p\ sequence (triangle), but each sequence shows a discontinuity, 
when the region of apparent intersection is well-resolved. (Figure 3 ofYoshida and 
Eriguchi, MNRAS [332]; used with permission.) 
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Figure 5: Frequencies of several axisymmetric modes along a sequence of rapidly ro- 
tating relativistic polytropes ofN= 1.0, in the Cowling approximation. On the hori- 
zontal axis, the angular velocity of each model is scaled to the angular velocity of the 
model at the mass-shedding limit. Lower-order modes are weakly affected by rapid 
rotation, while higher-order modes show apparent mode intersections. (Figure 10 of 
Font, Dimmelmeier, Gupta and Stergioulas, MNRAS [101].) 
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Quasi-radial modes in rotating relativistic stars have been studied by Hartle & 
Friedman [142] and by Datta et al. [22] in the slow-rotation approximation. Yoshida 
& Eriguchi [332] study quasi-radial modes of rapidly rotating stars in the relativistic 
Cowling approximation and find that apparent intersections between quasi-radial and 
other axisymmetric modes can appear near the mass-shedding Umit (see Figure 4). 
These apparent intersections are due to avoided crossings between mode sequences, 
which are also known to occur for axisymmetric modes of rotating Newtonian stars. 
Along a continuous sequence of computed mode frequencies an avoided crossing oc- 
curs when another sequence is encountered. In the region of the avoided crossing the 
eigenfunctions of the two modes become of mixed character. Away from the avoided 
crossing and along the continuous sequences of computed mode frequencies, the eigen- 
functions are exchanged. However, each "quasi-normal mode" is characterized by the 
shape of its eigenfunction and thus, the sequences of computed frequencies that belong 
to particular quasi-normal modes are discontinuous at avoided crossings (see Figure 
4 for more details). The discontinuities can be found in numerical calculations, when 
quasi-normal mode sequences are well resolved in the region of avoided crossings. 
Otherwise, quasi-normal mode sequences will appear as intersecting. 

Several axisymmetric modes have recently been computed for rapidly rotating rel- 
ativistic stars in the Cowhng approximation, using time-evolutions of the nonlinear 
hydrodynamical equations [101] (see [103] for a description of the 2D numerical evo- 
lution scheme). As in [332], Font et al. find that apparent mode intersections are 
common for various higher-order axisymmetric modes (see Figure 5). Axisymmetric 
inertial modes also appear in the numerical evolutions. 

The first fully relativistic frequencies of quasi-radial modes for rapidly rotating 
stars (without assuming the Cowling approximation) have been obtained recently, again 
through nonhnear time-evolutions [102] (see Section 4.2). 

• Going further The stabilization of a relativistic star, that is marginally stable 
to axisymmetric perturbations, by an external gravitational field, is discussed in 
[306]. 

3.4 Nonaxisymmetric Perturbations 
3.4.1 Nonrotating Limit 

Thorne, Campolattaro and Price, in a series of papers [310, 244, 305], initiated the com- 
putation of nonradial modes by formulating the problem in the Regge-Wheeler (RW) 
gauge [252] and numerically computing nonradial modes for a number of neutron star 
models. A variational method for obtaining eigenfrequencies and eigenfunctions has 
been constructed by Detweiler and Ipser [83]. Lindblom and Detweiler [199] explicitly 
reduced the system of equations to four first-order ordinary differential equations and 
obtained more accurate eigenfrequencies and damping times for a larger set of neutron 
star models. They later realized that their system of equations is sometimes singular in- 
side the star and obtained an improved set of equations which is free of this singularity 
[82]. 
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Chandrasekhar and Ferrari [59] expressed the nonradial pulsation problem in terms 
of a fifth-order system in a diagonal gauge, which is formally independent of fluid 
variables. They thus reformulate the problem in a way analogous to the scattering of 
gravitational waves off a black hole. Ipser and Price [160] show that in the RW gauge, 
nonradial pulsations can be described by a system of two second-order differential 
equations, which can also be independent of fluid variables. In addition, they find 
that the diagonal gauge of Chandrasekhar and Ferrari has a remaining gauge freedom 
which, when removed, also leads to a fourth-order system of equations [245]. 

In order to locate purely outgoing-wave modes, one has to be able to distinguish the 
outgoing-wave part from the ingoing-wave part at infinity. This is typically achieved 
using analytic approximations of the solution at infinity. 

W-modes pose a more challenging numerical problem because they are strongly 
damped and the techniques used for /— and p— modes fail to distinguish the outgoing- 
wave part. First accurate numerical solutions were obtained by Kokkotas and Schutz 
[178], followed by Leins, Nollert and Soffel [194]. Andersson, Kokkotas and Schutz 
[13], successfully combine a redefinition of variables with a complex-coordinate inte- 
gration method, obtaining highly accurate complex frequencies for w modes. In this 
method, the ingoing and outgoing solutions are separated by numerically calculating 
their analytic continuations to a place in the complex-coordinate plane, where they 
have comparable ampUtudes. Since this approach is purely numerical, it could prove to 
be suitable for the computation of quasi-normal modes in rotating stars, where analytic 
solutions at infinity are not available. 

The non-availabiUty of asymptotic solutions at infinity in the case of rotating stars is 
one of the major difficulties for computing outgoing modes in rapidly rotating relativis- 
tic stars. A method that may help to overcome this problem, at least to an acceptable 
approximation, has been found by Lindblom, Mendell and Ipser [200]. 

The authors obtain approximate near-zone boundary conditions for the outgoing 
modes that replace the outgoing-wave condition at infinity and that enable one to com- 
pute the eigenfrequencies with very satisfactory accuracy. First, the pulsation equations 
of polar modes in the Regge-Wheeler gauge are reformulated as a set of two second- 
order radial equations for two potentials - one corresponding to fluid perturbations and 
the other to the perturbations of the spacetime. The equation for the spacetime per- 
turbation reduces to a scalar wave equation at infinity and to Laplace's equation for 
zero-frequency solutions. From these, an approximate boundary condition for outgo- 
ing modes is constructed and imposed in the near zone of the star (in fact on its surface) 
instead of at infinity. For polytropic models, the near-zone boundary condition yields 
/-mode eigenfrequencies with real parts accurate to 0.01 % — 0. 1 % and imaginary parts 
with accuracy at the 10% — 20% level, for the most relativistic stars. If the near zone 
boundary condition can be applied to the oscillations of rapidly rotating stars, the re- 
sulting frequencies and damping times should have comparable accuracy. 
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3.4.2 Slow Rotation Approximation 

The slow rotation approximation is useful for obtaining a first estimate of the effect 
of rotation on the pulsations of relativistic stars. To lowest order in rotation, a polar 
/-mode of an initially nonrotating star couples to an axial Z ± 1 mode in the presence 
of rotation. Conversely, an axial /-mode couples to a polar / ± 1 mode as was first 
discussed by Chandrasekhar and Ferrari [59]. 

The equations of nonaxisymmetric perturbations in the slow -rotation limit are de- 
rived in a diagonal gauge, by Chandrasekhar and Ferrari [59], and in the Regge- 
Wheeler gauge, by Kojima [170, 172], where the complex frequencies a = aR + ioi 
for the / = m modes of various polytropes are computed. For counterrotating modes, 
both csr and a/ decrease, tending to zero, as the rotation rate increases (when a passes 
through zero, the star becomes unstable to the CFS-instabiUty). Extrapolating Or and 
O/ to higher rotation rates, Kojima finds a large discrepancy between the points where 
Or and O/ go through zero. This shows that the slow rotation formalism cannot ac- 
curately determine the onset of the CFS-instabiUty of polar modes in rapidly rotating 
neutron stars. 

In [171], it is shown that, for slowly rotating stars, the coupling between polar and 
axial modes affects the frequency of /- and p-vaoAes only to second order in rotation, 
so that, in the slow rotation approximation, to 0(fl), the coupUng can be neglected 
when computing frequencies. 

The linear perturbation equations in the slow-rotation approximation have recently 
been derived in a new gauge by Ruoff, Stavridis and Kokkotas [259]. Using the ADM 
formalism, a first-order hyperbolic evolution system is obtained, which is suitable for 
numerical integration without further manipulations (as was required in the Regge- 
Wheeler gauge). In this gauge (which is related to a gauge introduced for nonrotating 
stars in [27]) the symmetry between the polar and axial equations becomes directly 
apparent. 

The case of relativistic inertial modes is different, as these modes have both axial 
and polar parts at order 0{Q.) and the presence of continuous bands in the spectrum 
(at this order in the rotation rate) has lead to a series of detailed investigations of the 
properties of these modes (see [177] for a review). In a recent paper, Ruoff, Stavridis 
and Kokkotas [258] finally show that the inclusion of both polar and axial parts in 
the computation of relativistic r-modes, at order 0{D), allows for discrete modes to be 
computed, in agreement with post-Newtonian [212] and nonlinear, rapid-rotation [295] 
calculations. 

3.4.3 Post-Newtonian Approximation 

A step towards the solution of the perturbation equations in full general relativity has 
been taken by Cutler and Lindblom [75, 77, 196] who obtain frequencies for the / = 
m /-modes in rotating stars in the first post-Newtonian (1-PN) approximation. The 
perturbation equations are derived in the post-Newtonian formalism (see [133]), i.e. 
the equations are separated into equations of consistent order in 1/c. 
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Cutler and Lindblom show that in this scheme, the perturbation of the 1-PN cor- 
rection of the four- velocity of the fluid can be obtained analytically in terms of other 
variables; this is similar to the perturbation in the three-velocity in the Newtonian Ipser- 
Managan scheme. The perturbation in the 1-PN corrections are obtained by solving an 
eigenvalue problem, which consists of three second order equations, with the 1-PN 
correction to the eigenfrequency of a mode, Aco, as the eigenvalue. 

Cutler and Lindblom obtain a formula that yields Aco if one knows the 1-PN sta- 
tionary solution and the solution to the Newtonian perturbation equations. Thus, the 
frequency of a mode in the 1-PN approximation can be obtained without actually solv- 
ing the 1-PN perturbation equations numerically. The 1-PN code was checked in the 
nonrotating limit and it was found to reproduce the exact general relativistic frequen- 
cies for stars with M/R = 0.2, obeying an AT = 1 polytropic EOS, with an accuracy of 
3% -8%. 

Along a sequence of rotating stars, the frequency of a mode is commonly described 
by the ratio of the frequency of the mode in the comoving frame to the frequency of 
the mode in the nonrotating limit. For an A' = 1 polytrope and for M/R = 0.2, this 
frequency ratio is reduced by as much as 12% in the 1-PN approximation compared to 
its Newtonian counterpart (for the fundamental l = m modes) which is representative of 
the effect that general relativity has on the frequency of quasi-normal modes in rotating 
stars. 

3.4.4 Cowling Approximation 

In several situations, the frequency of pulsations in relativistic stars can be estimated 
even if one completely neglects the perturbation in the gravitational field, i.e. if one 
sets 5gab = in the perturbation equations 1224]. In this approximation, the pulsations 
are described only by the perturbation in the fluid variables and the scheme works quite 
well for /, p and r-modes [207]. A different version of the CowUng approximation, 
in which 5^,^ is kept nonzero in the perturbation equations, has been suggested to be 
more suitable for g-modes [96], since these modes could have large fluid velocities, 
even though the variation in the gravitational field is weak. 

Yoshida and Kojima [333], examine the accuracy of the relativistic Cowling ap- 
proximation in slowly rotating stars. The first-order correction to the frequency of a 
mode depends only on the eigenfrequency and eigenfunctions of the mode in the ab- 
sence of rotation and on the angular velocity of the star. The eigenfrequencies of /, pi 
and p2 modes for slowly rotating stars with M/R between 0.05 and 0.2 are computed 
(assuming polytropic EOSs with N =1 and A'^ = 1.5) and compared to their counter- 
parts in the slow- rotation approximation. 

For the Z = 2 /-mode, the relative error in the eigenfrequency because of the Cowl- 
ing approximation is 30% for weakly relativistic stars (M/R = 0.05) and about 15% 
for stars with M/R ~ 0.2 and the error decreases for higher /-modes. For the pi and 
P2 modes the relative error is similar in magnitude but it is smaller for less relativis- 
tic stars. Also, for p-modes, the Cowling approximation becomes more accurate for 
increasing radial mode number. 
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Figure 6: Eigenfrequencies (in the Cowling approximation) of f-modes along a M = 
I.SMq sequence of models, constructed with the WFF3-NV EOS. The vertical line 
corresponds to the frequency of rotation of the model at the mass-shedding limit of the 
sequence. (Figure 1 ofYoshida and Eriguchi, ApJ [331]; used with permission.) 



As an application, Yoshida and Eriguchi [330, 331] use the CowUng approximation 

to estimate the onset of the /-mode CFS instability in rapidly rotating relativistic stars 
and to compute frequencies of /-modes for several realistic equations of state (see 
Figure 6). 

3.5 Nonaxisymmetric Instabilities 
3.5.1 Introduction 

Rotating cold neutron stars, detected as pulsars, have a remarkably stable rotation pe- 
riod. But, at birth, or during accretion, rapidly rotating neutron stars can be subject to 
various nonaxisymmetric instabiUties, which will affect the evolution of their rotation 
rate. 

If a protoneutron star has a sufficiently high rotation rate (so that e.g. T /W > 
0.27 in the case of Maclaurin spheroids), it will be subject to a dynamical instability 
driven by hydrodynamics and gravity. Through the 1 = 2 mode, the instability will 
deform the star into a bar shape. This highly nonaxisymmetric configuration will emit 
strong gravitational waves with frequencies in the kHz regime. The development of the 
instability and the resulting waveform have been computed numerically in the context 
of Newtonian gravity by Houser et al. [152] and in full general relativity by Shibata et 
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al. [275] (see Section 4.1.3). 

At lower rotation rates, the star can become unstable to secular nonaxisymmetric in- 
stabilities, driven by gravitational radiation or viscosity. Gravitational radiation drives 
a nonaxisymmetric instability when a mode that is retrograde in a frame corotating 
with the star appears as prograde to a distant inertial observer, via the Chandrasekhar- 
Friedman-Schutz (CFS) mechanism [58, 115]: A mode that is retrograde in the corotat- 
ing frame has negative angular momentum, because the perturbed star has less angular 
momentum than the unperturbed one. If, for a distant observer, the mode is prograde, 
it removes positive angular momentum from the star and thus the angular momentum 
of the mode becomes increasingly negative. 

The instability evolves on a secular timescale, during which the star loses angular 
momentum via the emitted gravitational waves. When the star rotates slower than 
a critical value, the mode becomes stable and the instability proceeds on the longer 
timescale of the next unstable mode, unless it is suppressed by viscosity. 

Neglecting viscosity, the CFS -instability is generic in rotating stars for both polar 
and axial modes. For polar modes, the instability occurs only above some critical 
angular velocity, where the frequency of the mode goes through zero in the inertial 
frame. The critical angular velocity is smaller for increasing mode number I. Thus, 
there will always be a high enough mode number /, for which a slowly rotating star 
will be unstable. Many of the hybrid inertial modes (and in particular the relativistic 
r-mode) are generically unstable in all rotating stars, since the mode has zero frequency 
in the inertial frame when the star is nonrotating [4, 1 14]. 

The shear and bulk viscosity of neutron star matter is able to suppress the growth of 
the CFS-instability except when the star passes through a certain temperature window. 
In Newtonian gravity, it appears that the polar mode CFS-instability can occur only in 
nascent neutron stars that rotate close to the mass-shedding limit [158, 157, 154, 328, 
201], but the computation of neutral /-modes in full relativity [294, 297] shows that 
relativity enhances the instability, allowing it to occur in stars with smaller rotation 
rates than previously thought. 

• Going further. A numerical method for the analysis of the ergosphere instability 

in relativistic stars, which could be extended to nonaxisymmetric instabiUties of 
fluid modes, is presented by Yoshida and Eriguchi in [329]. 

3.5.2 CFS-InstabiUty of Polar Modes 

The existence of the CFS-instabiUty in rotating stars was first demonstrated by Chan- 
drasekhar [58] in the case of the 1 = 2 mode in uniformly rotating, uniform density 
Maclaurin spheroids. Friedman and Schutz [115] show that this instability also ap- 
pears in compressible stars and that all rotating self-gravitating perfect fluid configura- 
tions are generically unstable to the emission of gravitational waves. In addition, they 
find that a nonaxisymmetric mode becomes unstable when its frequency vanishes in 
the inertial frame. Thus, zero-frequency outgoing-modes in rotating stars are neutral 
(marginally stable). 
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In the Newtonian limit, neutral modes have been determined for several polytropic 
EOSs [153, 216, 156, 328]. The instability first sets in through / = m modes. Modes 
with larger I become unstable at lower rotation rates but viscosity limits the interesting 
ones to / < 5. For an A' = 1 polytrope, the critical values of T /W for the / = 3,4 
and 5 modes are 0.079, 0.058 and 0.045 respectively and these values become smaller 
for softer polytropes. The I = m = 2 "bar" mode has a critical T /W ratio of 0.14 
that is almost independent of the polytropic index. Since soft EOSs cannot produce 
models with high T /W values, the bar mode instability appears only for stiff Newtonian 
polytropes ofN< 0.808 [161, 284]. In addition, the viscosity-driven bar mode appears 
at the same critical T /W ratio as the bar mode driven by gravitational radiation [159] 
(we will see later that this is no longer true in general relativity). 

The post-Newtonian computation of neutral modes by Cutler and Lindblom [77, 
196] has shown that general relativity tends to strengthen the CFS-instability. Com- 
pared to their Newtonian counterparts, critical angular velocity ratios £lc/^o (where 
fio = (3Mo/4/?q)'/^ and Mq, Rq are the mass and radius of the nonrotating star in the 
sequence), are lowered by as much as 10% for stars obeying the = 1 polytropic 
EOS (for which the instabiUty occurs only for / = m > 3 modes in the post-Newtonian 
approximation). 

In full general relativity, neutral modes have been determined for polytropic EOSs 
of N > 1.0 by Stergioulas and Friedman [294, 297], using a new numerical scheme. 
The scheme completes the Eulerian formalism developed by Ipser & Lindblom in the 
Cowling approximation (where 5gab was neglected) [154], by finding an appropriate 
gauge in which the time-independent perturbation equations can be solved numerically 
for 8gab- The computation of neutral modes for polytropes of N = 1.0, 1.5 and 2.0 
shows that relativity significantly strengthens the instability. For the = 1 .0 polytrope, 
the critical angular velocity ratio Q.c/Q.k, where Q.k is the angular velocity at the mass- 
shedding limit at same central energy density, is reduced by as much as 15% for the 
most relativistic configuration (see Figure 7). A surprising result (which was not found 
in computations that used the post-Newtonian approximation), is that the I = m = 2 
bar mode is unstable even for relativistic polytropes of index = 1.0. The classical 
Newtonian result for the onset of the bar mode instability {Ncrit < 0.808) is replaced by 

Ncrit < 1.3, (54) 

in general relativity. It is evident that, in relativistic stars, the onset of the gravitational- 
radiation-driven bar mode does not coincide with the onset of the viscosity-driven bar 
mode, which occurs at larger T/W [38]. The computation of the onset of the CFS- 
instability in the relativistic Cowling approximation by Yoshida and Eriguchi [330] 
agrees qualitatively with the conclusions in [294, 297]. 

Morsink, Stergioulas and Blattning [228] extend the method presented in [297] to 
a wide range of realistic equations of state (which usually have a stiff high density re- 
gion, corresponding to polytropes of index A' = 0.5 — 0.7) and find that the / = m = 2 
bar mode becomes unstable for stars with gravitational mass as low as 1.0 — 1.2M0. 
For 1.4M0 neutron stars, the mode becomes unstable at 80% — 95% of the maximum 
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Figure 7: The I = m neutral f-mode sequences for EOS A. Shown are the ratio of rota- 
tional to gravitational energy T /W (upper panel) and the ratio of the critical angular 
velocity £lc to the angular velocity at the mass-shedding limit for uniform rotation 
(lower panel) as a function of gravitational mass. The solid curves are the neutral 
mode sequences for I = m = 2,3,4 and 5 (from top to bottom), while the dashed curve 
in the upper panel corresponds to the mass-shedding limit for uniform rotation. The 
I = m = 2 f-mode becomes CFS-unstable even at 85% of the mass-shedding limit, 
for 1.4 Mq models constructed with this EOS. (Figure 2 ofMorsink, Stergioulas and 
Blattning, ApJ [228].) 
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Figure 8: Eigenfrequencies (in the Cowling approximation) of the m — 2 mode as a 
function of the parameter P = for three different sequences of differentially ro- 

tating neutron stars (the A^^ = 0.0 line corresponding to uniform rotation). The filled 
dot indicates the neutral stability point of a uniformly rotating star computed in full 
General Relativity (Stergioulas and Friedman, 1998 [297]). Differential rotation shifts 
the neutral point to higher rotation rates. (Figure 1 ofYoshida, Rezzolla, Karino and 
Eriguchi inApJ[335]; used with permission.) 
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allowed rotation rate. For a wide range of equations of state, the / = m = 2 /-mode be- 
comes unstable at a ratio of rotational to gravitational energies T /W 0.08 for 1 AMq 
stars and T /W ~ 0.06 for maximum mass stars. This is to be contrasted with the New- 
tonian value ofT/W ~ 0.14. The empirical formula 



where M^ax is the maximum mass for a spherical star allowed by a given equation of 
state, gives the critical value of T /W for the bar /—mode instability, with an accuracy 
of 4% — 6%, for a wide range of realistic EOSs. 

In newly-born neutron stars the CFS-instability could develop while the back- 
ground equilibrium star is still differentially rotating. In that case, the critical value 
of T /W, required for the instability in the /-mode to set in, is larger than the corre- 
sponding value in the case of uniform rotation [335] (Figure 8). The mass-shedding 
limit for differentially rotating stars also appears at considerably larger T /W than the 
mass-shedding limit for uniform rotation. Thus, Yoshida et al. [335] suggest that 
differential rotation favours the instability, since the ratio (r/iy)criticai/(7'/I^)sheddmg 
decreases with increasing degree of differential rotation. 

3.5.3 CFS-InstabiUty of Axial Modes 

In nonrotating stars, axial fluid modes are degenerate at zero frequency but in rotating 
stars they have nonzero frequency and are called r-modes in the Newtonian limit [241, 
261]. To order their frequency in the inertial frame is 



while the radial eigenfunction of the perturbation in the velocity can be determined at 
order [173]. According to (56), r-modes with m > are prograde (co,- < 0) with 
respect to a distant observer but retrograde (co, = to, + m£l > 0) in the comoving frame 
for all values of the angular velocity. Thus, r-modes in relativistic stars are generically 
unstable to the emission of gravitational waves via the CFS-instability, which was first 
discovered by Andersson [4], in the case of slowly rotating, relativistic stars. This 
result was proved rigorously by Friedman and Morsink [114], who showed that the 
canonical energy of the modes is negative. 

Two independent computations in the Newtonian Cowling approximation [206, 14] 
showed that the usual shear and bulk viscosity assumed to exist for neutron star matter 
is not able to damp the r-mode instability, even in slowly rotating stars. In a temperature 
window of 10^ K <T < 10^*^ K, the growth time of the I = m = 2 mode becomes 
shorter than the shear or bulk viscosity damping time at a critical rotation rate that is 
roughly one tenth the maximum allowed angular velocity of uniformly rotating stars. 
The gravitational radiation is dominated by the mass current quadrupole term. These 
results suggested that a rapidly rotating proto-neutron star will spin down to Crab-Uke 
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Figure 9: The r-mode instability window for a strange star of M = 1.4M© and R = 
10 km (solid line). Dashed curves show the corresponding instability windows for 
normal npe fluid and neutron stars with a crust. The instability window is compared 
to i) the inferred spin-periods for accreting stars in Low-mass X-ray binaries [shaded 
box], and ii) the fastest known millisecond pulsars (for which observational upper 
limits on the temperature are available) [horizontal lines]. (Figure 1 of Andersson, 
Jones and Kokkotas [9]; used with permission.) 

rotation rates within one year of its birth, because of the r-mode instability. Due to 
uncertainties in the actual viscosity damping times and because of other dissipative 
mechanisms, this scenario also is consistent with somewhat higher initial spins, like 
the suggested initial spin of a several milliseconds of the X-ray pulsar in the supernova 
remnant N157B [222]. Millisecond pulsars with periods less than a few milliseconds 
can then only form after the accretion-induced spin-up of old pulsars and not in the 
accretion-induced collapse of a white dwarf. 

The precise limit on the angular velocity of newly-born neutron stars will depend 
on several factors, such as the strength of the bulk viscosity, the cooling process, su- 
perfluidity, presence of hyperons, influence of a solid crust etc. In the uniform density 
approximation, the r-mode instability can be studied analytically to 0{Qr') in the an- 
gular velocity of the star [179]. A study on the issue of detectability of gravitational 
waves from the r-mode instability was presented in [237] (see section 3.5.5), while 
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Figure 10: Relativistic r-mode frequencies for a range of the compactness ratio M/R. 
The coupling of polar and axial terms, even in the order 0{D.) slow-rotation approx- 
imation has a dramatic impact on the continuous frequency bands (shaded areas), 
allowing the r-mode to exist even in highly compact stars. The Newtonian value of the 
r-mode frequency is plotted as a dashed-dotted line. (Figure 3 ofRuoff, Stavridis and 
Kokkotas [258]; used with permission.) 

Andersson, Kokkotas and Stergioulas [15] and Bildsten [35] proposed that the r-mode 
instability is limiting the spin of millisecond pulsars spun-up in LMXBs and it could 
even set the minimum observed spin period of ^ 1.5 ms (see [10]). This scenario is 
also compatible with observational data, if one considers strange stars instead of neu- 
tron stars [9] (see Figure 9). 

Since the discovery of the r-mode instability, a large number of authors have stud- 
ied in more detail the development of the instability and its astrophysical consequences. 
Unlike in the case of the /-mode instabiUty, many different aspects and interactions 
have been considered. This intense focus on the detailed physics has been very fruit- 
ful and we now have a much more complete understanding of the various physical 
processes that are associated with pulsations in rapidly rotating relativistic stars. The 
latest understanding of the r-mode instability is that it may not be a very promising 
gravitational wave source (as originally thought), but the important astrophysical con- 
sequences, such as the limits of the spin of young and of recycled neutron stars are 
still considered plausible. The most crucial factors affecting the instability are mag- 
netic fields [288, 256, 254, 255], possible hyperon bulk viscosity [163, 205, 137] and 
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Figure 1 1 : Projected trajectories of several fiducial fluid elements ( as seen in the coro- 
tating frame) for an I = m = 2 Newtonian r-mode. All of the fluid elements are initially 
positioned on the — Q meridian, at different latitudes (indicated with stars). Blue 
dots indicate the position of the fluid elements after each full oscillation period. The r- 
mode induces a kinematical, differential drift. (Figure 2c of Rezzolla, Lamb, Markovic 
and Shapiro, PRD [254]; used with permission.) 
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Figure 12: Evolution of the axial velocity in the equatorial plane for a relativistic r- 
mode in a rapidly rotating N = l.O polytrope (in the Cowling approximation). Since 
the initial data used to excite the mode are not exact, the evolution is a superposition 
of ( mainly) the I = m = 2 r-mode and several inertial modes. The amplitude of the 
oscillation decreases due to numerical (finite-differencing) viscosity of the code. A 
beating between the I = m = 2 r-mode and another inertial mode can also be seen. 
(Figure 2 of Stergioulas and Font, PRL [295].) 



nonlinear saturation [295, 208, 209, 19]. The question of the possible existence of a 
continuous spectrum has also been discussed by several authors, but the most recent 
analysis suggests that higher-order rotational effects still allow for discrete r-modes in 
relativistic stars [334, 258] (see Figure 10). 

Magnetic fields can affect the r-mode instability, as the r-mode velocity field cre- 
ates differential rotation, which is both kinematical and due to gravitational radiation 
reaction (see Figure 11). Under differential rotation, an initially weak poloidal mag- 
netic field is wound-up, creating a strong toroidal field, which causes the r-mode am- 
plitude to saturate. If neutron stars have hyperons in their cores, the associated bulk 
viscosity is so strong that it could completely prevent the growth of the r-mode in- 
stabiUty. However, hyperons are predicted only by certain equations of state and the 
relativistic mean field theory is not universally accepted. Thus, our ignorance of the 
true equation of state still leaves a lot of room for the r-mode instabihty to be consid- 
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ered viable. 

The gravitational-wave detection from r-modes depends crucially on the nonhnear 
saturation amplitude. A first study by Stergioulas & Font [295] suggests that r-modes 
can exist at large amplitudes of order unity for dozens of rotational periods in rapidly 
rotating relativistic stars (Figure 12). The study used 3D relativistic hydrodynamical 
evolutions in the Cowling approximation. This result was confirmed by Newtonian 
3D simulations of nonlinear r- modes by Lindblom, Tohline & Vallisneri [205, 208]. 
Lindblom et al. went further, using an accelerated radiation reaction force, to artificially 
grow the r-mode amplitude on a hydrodynamical (instead of the secular) timescale. At 
the end of the simulations, the r-mode grew so large that large shock waves appeared 
on the surface of the star, while the amplitude of the mode subsequently collapsed. 
Lindblom et al. suggested that shock heating may be the mechanism that saturates the 
r-modes at a dimensionless amplitude of a ~ 3. 

More recent studies of nonlinear couplings between the r-mode and higher-order 
inertial modes [19] and new 3D nonlinear Newtonian simulations [132] seem to sug- 
gest a different picture. The r-mode could be saturated due to mode couplings or due to 
a hydrodynamical instability at amplitudes much smaller than the amplitude at which 
shock waves appeared in the simulations by Lindblom et al. Such a low amplitude, on 
the other hand, modifies the properties of the r-mode instability as a gravitational wave 
source, but is not necessarily bad news for gravitational wave detection, as a lower 
spin-down rate also implies a higher event rate for the r-mode instability in LMXBs in 
our own Galaxy [9, 151]. The 3D simulations need to achieve significantly higher reso- 
lutions before definite conclusions can be reached, while the Arras et al. work could be 
extended to rapidly rotating relativistic stars (in which case the mode frequencies and 
eigenfunctions could change significantly, compared to the slowly rotating Newtonian 
case, which could affect the nonlinear coupling coefficients). Spectral methods can 
be used for achieving high accuracy in mode calculations and first results have been 
obtained by Villain and Bonazzolla [318] for inertial modes of slowly rotating stars in 
the relativistic Cowling approximation. 

For a more extensive coverage of the numerous articles on the r-mode instability 
that appeared in recent years, the reader is referred to several excellent recent review 
articles [12, 113, 198, 177, 5]. 

• Going further. If rotating stars with very high compactness exist, then w-modes 

can also become unstable, as was recently found by Kokkotas, Ruoff and Ander- 
sson [180]. The possible astrophysical imphcations are still under investigation. 

3.5.4 Effect of Viscosity on CFS-InstabUity 

In the previous sections, we have discussed the growth of the CFS-instability driven 
by gravitational radiation in an otherwise nondissipative star The effect of neutron 
star matter viscosity on the dynamical evolution of nonaxisymmetric perturbations can 
be considered separately, when the timescale of the viscosity is much longer than the 
oscillation timescale. If Xgr is the computed growth rate of the instabiUty in the absence 
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of viscosity and Xj, Xh are the timescales of shear and bulk viscosity, then the total 
timescale of the perturbation is 

l = -!- + i + l. (57) 

X Tor Xj Tb 

Since Xgr < and Zh, > 0, a mode will grow only if Xgr is shorter than the viscous 
timescales, so that 1/t < 0. 

In normal neutron star matter, shear viscosity is dominated by neutron-neutron scat- 
tering with a temperature dependence of T^^ |98| and computations in the Newtonian 
limit and post-Newtonian approximation show that the CFS -instability is suppressed 
for r < 10*^ K - 10^ K [158, 157, 328, 196]. If neutrons become a superfluid be- 
low a transition temperature Tg, then mutual friction, which is caused by the scattering 
of electrons off the cores of neutron vortices could significantly suppress the /-mode 
instability for T < [203] but the r-mode instability remains unaffected [204]. The 
superfluid transition temperature depends on the theoretical model for superfluidity and 
hes in the range 10** K - 6 x 10*^ K [239]. 

In a pulsating fluid that undergoes compression and expansion, the weak interac- 
tion requires a relatively long time to re-establish equilibrium. This creates a phase 
lag between density and pressure perturbations, which results in a large bulk viscosity 
[264]. The bulk viscosity due to this effect can suppress the CFS-instabiUty only for 
temperatures for which matter has become transparent to neutrinos [188, 40]. It has 
been proposed that for T > 5 x lO^K, matter will be opaque to neutrinos and the neu- 
trino phase space could be blocked ([188] see also [40]). In this case, bulk viscosity 
will be too weak to suppress the instabiUty, but a more detailed study is needed. 

In the neutrino transparent regime, the effect of bulk viscosity on the instability 
depends crucially on the proton fraction Xp. If Xp is lower than a critical value (~ 5), 
only modified URCA processes are allowed and bulk viscosity limits, but does not 
suppress completely, the instabihty [158, 157, 328]. For most modern EOSs, however, 
the proton fraction is larger than ~ ^ at sufficiently high densities [191], allowing 
direct URCA processes to take place. In this case, depending on the EOS and the 
central density of the star, the bulk viscosity could almost completely suppress the 
CFS-instabiUty in the neutrino transparent regime [339]. At high temperatures, T > 
5 X 10^ K, even if the star is opaque to neutrinos, the direct URCA cooling timescale 
to r 5 X 10^ K could be shorter than the the growth timescale of the CFS instability. 

3.5.5 Gravitational Radiation from CFS-Instability 

Conservation of angular momentum and the inferred initial period (assuming magnetic 
braking) of a few milliseconds for the X-ray pulsar in the supernova remnant N157B 
[222], suggests that a fraction of neutron stars may be born with very large rotational 
energies. The /-mode bar CFS-instabiUty thus appears as a promising source for the 
planned gravitational wave detectors [188]. It could also play a role in the rotational 
evolution of merged binary neutron stars, if the post-merger angular momentum ex- 
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ceeds the maximum allowed to form a Kerr black hole [28] or if differential rotation 
temporarily stabilizes the merged object. 

Lai and Shapiro [188] have studied the development of the /-mode instability us- 
ing Newtonian ellipsoidal models [186, 187]. They consider the case when a rapidly 
rotating neutron star is created in a core collapse. After a brief dynamical phase, the 
protoneutron star becomes secularly unstable. The instability deforms the star into a 
nonaxisymmetric configuration via the / = 2 bar mode. Since the star loses angular 
momentum via the emission of gravitational waves, it spins-down until it becomes sec- 
ularly stable. The frequency of the waves sweeps downward from a few hundred Hz to 
zero, passing through LIGO's ideal sensitivity band. A rough estimate of the wave am- 
plitude shows that, at lOOHz, the gravitational waves from the CFS-instability could 
be detected out to the distance of 140Mpc by the advanced LIGO detector This result is 
very promising, especially since for relativistic stars the instability will be stronger than 
the Newtonian estimate [297]. Whether /--modes should also be considered a promis- 
ing gravitational- wave source depends crucially on their nonlinear saturation amplitude 
(see Section 3.5.3). 

• Going further. The possible ways for neutron stars to emit gravitational waves 
and their detectability are reviewed in [43, 44, 1 17, 97, 307, 267, 78]. 

3.5.6 Viscosity-Driven Instability 

A different type of nonaxisymmetric instability in rotating stars is the instability driven 
by viscosity, which breaks the circulation of the fluid [257, 161]. The instability is 
suppressed by gravitational radiation, so it cannot act in the temperature window in 
which the CFS-instability is active. The instability sets in when the frequency of an 
/ = —m mode goes through zero in the rotating frame. In contrast to the CFS-instability, 
the viscosity-driven instability is not generic in rotating stars. The m = 2 mode becomes 
unstable at a high rotation rate for very stiff stars and higher m-modes become unstable 
at larger rotation rates. 

In Newtonian polytropes, the instability occurs only for stiff polytropes of index 
< 0.808 [161, 284]. For relativistic models, the situation for the instability becomes 
worse, since relativistic effects tend to suppress the viscosity-driven instability (while 
the CFS-instability becomes stronger). According to recent results by Bonazzola et al. 
[38], for the most relativistic stars, the viscosity-driven bar mode can become unstable 
only UN < 0.55. For 1.4M0 stars, the instabihty is present forA'^ < 0.67. 

These results are based on an approximate computation of the instability in which 
one perturbs an axisymmetric and stationary configuration and studies its evolution by 
constructing a series of triaxial quasi-equilibrium configurations. During the evolu- 
tion only the dominant nonaxisymmetric terms are taken into account. The method 
presented in [38] is an improvement (taking into account nonaxisymmetric terms of 
higher order) of an earlier method by the same authors [40]. Although the method 
is approximate, its results indicate that the viscosity-driven instability is likely to be 
absent in most relativistic stars, unless the EOS turns out to be unexpectedly stiff. 
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An investigation of the viscosity-driven bar mode instability, using incompressible, 
uniformly rotating triaxial ellipsoids in the post-Newtonian approximation, by Shapiro 
and Zane [270], finds that the relativistic effects increase the critical T /W ratio for the 
onset of the instability significantly. More recently, new post-Newtonian [85] and fully 
relativistic calculations for uniform-density stars [125] show that the viscosity-driven 
instability is not as strongly suppressed by relativistic effects as suggested in [270]. 
The most promising case for the onset of the viscosity-driven instability (in terms of 
the critical rotation rate) would be rapidly rotating strange stars [126], but the instability 
can only appear if its growth rate is larger than the damping rate due to the emission of 
gravitational radiation - a corresponding detailed comparison is still missing. 

4 Rotating Stars in Numerical Relativity 

Recently, the dynamical evolution of rapidly rotating stars has become possible in nu- 
merical relativity. In the framework of the 3+1 split of the Einstein equations [285] a 
stationary axisymmetric star can be described by a metric of the standard form 

ds^ = - (a^ - %^')dt^ + 2^idx'dt + ytjdx'dx^, (58) 

where a is the lapse function, P' is the shift three-vector and ytj is the spatial three- 
metric, with J = 1 ... 3. The spacetime has the following properties: 

• The metric function co in (5) describing the dragging of inertial frames by rotation 
is related to the shift vector through ^'^ = — co. This shift vector satisfies the 
minimal distortion shift condition. 

• The metric satisfies the maximal slicing condition, while the lapse function is 
related to the metric function v in (5) through a = e^. 

• The quasi-isotropic coordinates are suitable for numerical evolution, while the 
radial-gauge coordinates [25] are not suitable for non-spherical sources (see [46] 
for details). 

• The zero-angular momentum observers (ZAMOs) are the Eulerian observers, 
whose worldlines are normal to the ?=const. hypersurfaces. 

• Uniformly rotating stars have n=const. in the coordinate frame. This can be 
shown by requiring a vanishing rate of shear. 

• Normal modes of pulsation are discrete in the coordinate frame and their fre- 
quencies can be obtained by Fourier transforms (with respect to coordinate time 
t) of evolved variables at a fixed coordinate location [103]. 

Crucial ingredients for the successful long-term evolutions of rotating stars in nu- 
merical relativity are the conformal ADM schemes for the spacetime evolution (see 
[233, 278, 32, 2]) and hydrodynamical schemes that have been shown to preserve well 
the sharp rotational profile at the surface of the star [103, 295, 102]. 
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Figure 13: Time-evolution of the rotational velocity profile for a stationary, rapidly 
rotating relativistic star (in the Cowling approximation), using the 3rd order PPM 
scheme and a 116^ grid. The initial rotational profile is preserved to a high degree 
of accuracy, even after 20 rotational periods. (Figure 1 of Stergioulas and Font, PRL 
[295].) 

4.1 Numerical Evolution of Equilibrium Models 
4.1.1 Stable Equilibrium 

The long-term stable evolution of rotating relativistic stars in 3-D simulations has be- 
come possible through the use of High-Resolution-Shock-Capturing (HRSC) methods 
(see [100] for a review). Stergioulas and Font [295] evolve rotating relativistic stars 
near the mass-shedding limit for dozens of rotational periods (evolving only the equa- 
tions of hydrodynamics) (see Figure 13), while accurately preserving the rotational 
profile, using the 3rd order PPM method [63]. This method was shown to be superior 
to other, commonly used methods, in 2-D evolutions of rotating relativistic stars [103]. 

Fully coupled hydrodynamical and spacetime evolutions in 3-D have been obtained 
by Shibata [272] and by Font et al. [102]. In [272] the evolution of approximate 
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Figure 14; Simulation of a stationary, rapidly rotating neutron star model in full Gen- 
eral Relativity, for 3 rotational periods ( shown are iso-density contours, in dimension- 
less units - the figure only shows the initial configuration). The stationary shape is 
well-preserved at a resolution of 129^. Simulation by Font, Goodale, Iyer, Miller, Rez- 
zolla, Seidel, Stergioulas, Suen and Tobias. Visualization by W. Benger and L. Rezzolla 
at the Albert-Einstein-Institute, Golm [230]. 
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(conformally flat) initial data, is presented for about two rotational periods, and in [102] 
the simulations extend to several full rotational periods, using numerically exact initial 
data and a monotonized central difference (MC) slope limiter [317]. The MC slope 
hmiter is somewhat less accurate in preserving the rotational profile of equilibrium 
stars, than the 3rd order PPM method, but, on the other hand, it is easier to implement 
in a numerical code. 

New evolutions of uniformly and differentially rotating stars in 3-D, using different 
gauges and coordinate systems, are presented in [90], while new 2-D evolutions are 
presented in [271]. 

4.1.2 Instability to Collapse 

Shibata, Baumgarte and Shapiro [276] study the stability of supramassive neutron stars 
rotating at the mass-shedding limit, for a F = 2 polytropic EOS. Their 3-D simulations 
in full general relativity show that stars on the mass-shedding sequence with central 
energy density somewhat larger than that of the maximum mass model, are dynami- 
cally unstable to collapse. Thus, the dynamical instability of rotating neutron stars to 
axisyrmnetric perturbations is close to the corresponding secular instabihty. The initial 
data for these simulations are approximate, conformally flat axisymmetric solutions, 
but their properties are not very different from exact axisymmetric solutions even near 
the mass-shedding limit [71]. It should be noted that the approximate minimal distor- 
tion (AMD) shift condition does not prove useful in the numerical evolution, once a 
horizon forms. Instead, modified shift conditions are used in [276]. In the above simu- 
lations, no massive disk around the black hole is formed, as the equatorial radius of the 
initial model is inside the radius which becomes the ISCO of the final black hole. This 
could change if a different EOS is chosen. 

4.1.3 Dynamical Bar-Mode Instability 

Shibata, Baumgarte and Shapiro [275] study the dynamical bar-mode instabihty in 
differentially rotating neutron stars, in fuUy relativistic 3-D simulations. They find 
that stars become unstable when rotating faster than a critical value of P = T jW ~ 
0.24 — 0.25. This is only somewhat smaller than the Newtonian value of P ~ 0.27. 
Models with rotation only somewhat above critical become differentially rotating el- 
lipsoids, while models with P much larger than critical also form spiral arms, leading to 
mass ejection (see Figures 15, 16 and 17). In any case, the differentially rotating ellip- 
soids formed during the bar-mode instability have P > 0.2 indicating that they will be 
secularly unstable to bar-mode formation, driven by gravitational radiation or viscosity. 
The decrease of the critical value of P for dynamical bar formation due to relativistic 
effects has been confirmed by post-Newtonian simulations [260]. 
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Figure 15: Density contours and velocity flow for a neutron star model that has de- 
veloped spiral arms, due to the dynamical bar-mode instability. The computation was 
done in full General Relativity. (Figure 4 ofShibata, Baumgarte, and Shapiro in ApJ 
[275]; used with permission). 
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Figure 16: Simulation of the development of the dynamical bar-mode instability in a 
rapidly rotating relativistic star Spiral arms form within a few rotational periods. The 
different colors correspond to different values of the density, while the computation 
was done in full General Relativity. Movie produced at the University of Illinois by 
Thomas W. Baumgarte, Stuart L. Shapiro and Masaru Shibata, with the assistance of 
the Illinois Undergraduate Research Team [29]; used with permission. 
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Figure 17: Gravitational wave emission during the development of the dynamical bar- 
mode instability in a rapidly rotating relativistic star The gravitational wave amplitude 
in a plane containing the rotation axis is shown. At large distances, the waves assume 
a quadrupole-like angular dependence. Movie produced at the University of Illinois by 
Thomas W. Baumgarte, Stuart L. Shapiro and Masaru Shibata, with the assistance of 
the Illinois Undergraduate Research Team [29]; used with permission. 
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4.2 Pulsations of Rotating Stars 



Pulsations of rotating relativistic stars are traditionally studied (when possible) as a 
time-independent, linear eigenvalue problem, but recent advances in numerical rela- 
tivity also allow the study of such pulsations via numerical time-evolutions. The first 
quasi-radial mode frequencies of rapidly rotating stars in full general relativity have 
been recently obtained in [102], something that has not been achieved yet with the lin- 
ear perturbation theory. The fundamental quasi-radial mode in full general relativity 
has a similar rotational dependence as in the relativistic Cowling approximation and 
an empirical relation between the full GR computation and the Cowling approxima- 
tion can be constructed (Figure 18). For higher-order modes, apparent intersections of 
mode sequences near the mass-shedding limit do not allow for such empirical relations 
to be constructed. 

In the relativistic Cowling approximation, 2-D time-evolutions have yielded fre- 
quencies for the / = to / = 3 axisymmetric modes of rapidly rotating relativistic 
polytropes of A'^ = 1.0 [101]. The higher-order overtones of these modes show char- 
acteristic apparent crossings near mass-shedding (as was observed for the quasi-radial 
modes in [332]). 

Numerical relativity has also enabled the first study of nonlinear r-modes in rapidly 
rotating relativistic stars (in the Cowling approximation) by Stergioulas and Font [295]. 
For several dozen dynamical timescales, the study shows that nonlinear r-modes with 
amplitudes of order unity can exist in a star rotating near mass-shedding. However, on 
longer timescales, nonlinear effects may limit the r-mode amplitude to smaller values 
(see Section 3.5.3). 

4.3 Rotating Core Collapse 
4.3.1 Collapse to a Rotating Black Hole 

Black hole formation in relativistic core collapse was first studied in axisymmetry by 
Nakamura [231, 232], using the (2h-1)h-1 formalism [215]. The outcome of the simula- 
tion depends on the rotational parameter 

q = J/M^. (59) 

A rotating black hole is formed, only \f q < 1, indicating that cosmic censorship holds. 
Stark and Piran [29 1, 243] use the 2>+ 1 formalism and the radial gauge of Bardeen-Piran 
[25] to study black hole formation and gravitational wave emission in axisynnmetry. In 
this gauge, two metric functions used in determining gee and gi/ts/ can be chosen such 
that they tend directly at large radii to h+ and (the even and odd transverse traceless 
amplitudes of the gravitational waves, with \ /r fall-off at large radii) ^ In this way, the 
gravitational waveform is obtained at large radii directly in the numerical evolution. It 
is also easy to compute the gravitational energy emitted, as a simple integral over a 
sphere far from the source: AE ^ r^ / dt{h\ ,. + h^ ,.)- Using polar slicing, black hole 

defined in [291] has the opposite sign as that commonly used, e.g. in [309] 
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Figure 18: The first fiilly-relativistic, quasi-radial pulsation frequencies for a sequence 
of rapidly rotating stars (solid lines). The frequencies of the fundamental mode F 
(filled squares) and of the first overtone H\ (filled circles) are obtained through coupled 
hydrodynamical and spacetime evolutions. The corresponding frequencies obtained 
from computations in the relativistic Cowling approximation (in [101]) are shown as 
dashed lines. (Figure 16 of Font, Goodale, Iyer, Miller, Rezzolla, Seidel, Stergioulas, 
Suen and Tobias, PRD [102].) 
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formation appears as a region of exponentially small lapse, when ^ < 0(1). The ini- 
tial data consists of a nonrotating, pressure deficient, TOV solution, to which angular 
momentum is added by hand. The obtained waveform is nearly independent of the 
details of the collapse: it consists of a broad initial peak (since the star adjusts its initial 
spherical shape to a flattened shape, more consistent with the prescribed angular mo- 
mentum), the main emission (during the formation of the black hole) and an oscillatory 
tail, corresponding to oscillations of the formed black hole spacetime. The energy of 
the emitted gravitational waves during the axisymmetric core coUapse is found not to 
exceed 7 x IO^^Mqc^ (to which the broad initial peak has a negligible contribution). 
The emitted energy scales as q"^, while the energy in the even mode exceeds by at least 
an order of magnitude that in the odd mode. 

More recently, Shibata [274] carried out axisymmetric simulations of rotating stel- 
lar collapse in full general relativity, using a Cartesian grid, in which axisymmetry 
is imposed by suitable boundary conditions. The details of the formalism (numerical 
evolution scheme and gauge) are given in [273]. It is found that rapid rotation can 
prevent prompt black hole formation. When q = 0(1), a prompt collapse to a black 
hole is prevented even for a rest mass that is 70-80% larger than the maximum allowed 
mass of spherical stars and this depends weakly on the rotational profile of the ini- 
tial configuration. The final configuration is supported against collapse by the induced 
differential rotation. In these axisymmetric simulations, shock formation for q < 0.5 
does not result in a significant heating of the core: shocks are formed at a spheroidal 
shell around the high density core. In contrast, when the initial configuration is rapidly 
rotating (^ = 0(1 )), shocks are formed in a highly nonspherical manner near high den- 
sity regions and the resultant shock heating contributes in preventing prompt collapse 
to a black hole. A qualitative analysis in [274] suggests that a disk can form around 
a black hole, during core collapse, provided the progenitor is nearly rigidly rotating 
and q = 0(1) for a stiff progenitor EOS. On the other hand, ^ << 1 still allows for 
a disk formation, if the progenitor EOS is soft. At present, it is not clear how much 
the above conclusions depend on the restriction to axisymmetry or on other assump- 
tions - 3 -dimensional simulations of the core coUapse of such initiaUy axisymmetric 
configurations have still to be performed. 

A new numerical code for axisymmetric gravitational collapse in the (2-h1)-h1 for- 
mahsm is presented in [61]. 

4.3.2 Formation of Rotating Neutron Stars 

First attempts in studying the formation of rotating neutron stars in axisymmetric col- 
lapse were initiated by Evans [93, 94]. Recently, Dimmelmeier, Font & Miiller [89, 86] 
have successfully obtained detailed simulations of neutron star formation in rotating 
collapse. In the numerical scheme, HRSC methods are employed for the hydrodynam- 
ical evolution, while for the spacetime evolution the conformal flatness approximation 
[326] is used. Surprisingly, the gravitational waves obtained during the neutron star 
formation in rotating core coUapse are weaker in general relativity than in Newtonian 
simulations. The reason for this result is that relativistic rotating cores bounce at larger 
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central densities than in the Newtonian limit (for the same initial conditions). The 
gravitational waves are computed from the time-derivatives of the quadrupole moment, 
which involves the volume integration of pr^. As the density profile of the formed neu- 
tron star is more centrally condensed than in the Newtonian case, the corresponding 
gravitational waves turn out to be weaker. Details of the numerical methods and of the 
gravitational wave extraction used in the above studies can be found in [87, 88]. 

New, fully relativistic axisymmetric simulations, with coupled hydrodynamicaland 
spacetime evolution in the light-cone approach, have been obtained by Siebel et al. 
[283, 282] One of the advantages of the light-cone approach is that gravitational waves 
can be extracted accurately at null infinity, without spurious contamination by bound- 
ary conditions. The code by Siebel et al. combines the light-cone approach for the 
spacetime evolution with HRSC methods for the hydrodynamical evolution. In [282] 
it is found that gravitational waves are extracted more accurately using the Bondi news 
function, than by a quadrupole formula on the null cone. 

A new 2-D code for axisymmetric core collapse, also using HRSC methods has 
recently been introduced in [271]. 
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